lars/FYS-STK4155
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge pull request #4 from cforssen/bayesian

Browse Source 
More text in the Bayesian chapter.
This commit is contained in:
Morten Hjorth-Jensen
2018-03-13 19:29:33 +01:00
committed by GitHub
parent 6a1d62ceae 82795f8761
commit e7bf07a13b
20 changed files with 2784 additions and 60 deletions
Download Patch File Download Diff File Expand all files Collapse all files
doc/pub/Bayesian/html/._Bayesian-bs000.html
+52 -12
View File
@@ -48,7 +48,35 @@ Automatically generated HTML file from DocOnce source
("Bayes' theorem", 2, None, '___sec5'),
("The friends of Bayes' theorem", 2, None, '___sec6'),
('Inference With Parametric Models', 2, None, '___sec7'),
('Illustrative examples with python code', 2, None, '___sec8')]}
('Illustrative examples with python code', 2, None, '___sec8'),
('Example: Is this a fair coin?', 2, None, '___sec9'),
('A few words on different priors', 2, None, '___sec10'),
('Bayesian parameter estimation (single parameter)',
2,
None,
'___sec11'),
('Example: Measured flux from a star', 3, None, '___sec12'),
('Simple Photon Counts: Frequentist Approach',
3,
None,
'___sec13'),
('Simple Photon Counts: Bayesian Approach', 3, None, '___sec14'),
('A note about priors', 3, None, '___sec15'),
('Simple Photon Counts: Bayesian approach in practice',
3,
None,
'___sec16'),
('Best estimates and confidence intervals', 3, None, '___sec17'),
('Simple Photon Counts: Best estimates and confidence intervals',
3,
None,
'___sec18'),
('Bayesian parameter estimation (multiple parameters, '
'covariance)',
2,
None,
'___sec19'),
('Bayesian model selection', 2, None, '___sec20')]}
end of tocinfo -->
<body>
@@ -86,15 +114,27 @@ MathJax.Hub.Config({
<li class="dropdown">
<a href="#" class="dropdown-toggle" data-toggle="dropdown">Contents <b class="caret"></b></a>
<ul class="dropdown-menu">
<!-- navigation toc: --> <li><a href="._Bayesian-bs001.html#___sec0" style="font-size: 80%;">What is Bayesian Statistics</a></li>
<!-- navigation toc: --> <li><a href="._Bayesian-bs002.html#___sec1" style="font-size: 80%;">Inference</a></li>
<!-- navigation toc: --> <li><a href="._Bayesian-bs003.html#___sec2" style="font-size: 80%;">Statistical Inference</a></li>
<!-- navigation toc: --> <li><a href="._Bayesian-bs004.html#___sec3" style="font-size: 80%;">Some history</a></li>
<!-- navigation toc: --> <li><a href="._Bayesian-bs010.html#___sec4" style="font-size: 80%;">The Bayesian recipe</a></li>
<!-- navigation toc: --> <li><a href="._Bayesian-bs011.html#___sec5" style="font-size: 80%;">Bayes' theorem</a></li>
<!-- navigation toc: --> <li><a href="._Bayesian-bs013.html#___sec6" style="font-size: 80%;">The friends of Bayes' theorem</a></li>
<!-- navigation toc: --> <li><a href="._Bayesian-bs014.html#___sec7" style="font-size: 80%;">Inference With Parametric Models</a></li>
<!-- navigation toc: --> <li><a href="._Bayesian-bs016.html#___sec8" style="font-size: 80%;">Illustrative examples with python code</a></li>
<!-- navigation toc: --> <li><a href="._Bayesian-bs001.html#___sec0" style="font-size: 80%;"><b>What is Bayesian Statistics</b></a></li>
<!-- navigation toc: --> <li><a href="._Bayesian-bs002.html#___sec1" style="font-size: 80%;"><b>Inference</b></a></li>
<!-- navigation toc: --> <li><a href="._Bayesian-bs003.html#___sec2" style="font-size: 80%;"><b>Statistical Inference</b></a></li>
<!-- navigation toc: --> <li><a href="._Bayesian-bs004.html#___sec3" style="font-size: 80%;"><b>Some history</b></a></li>
<!-- navigation toc: --> <li><a href="._Bayesian-bs010.html#___sec4" style="font-size: 80%;"><b>The Bayesian recipe</b></a></li>
<!-- navigation toc: --> <li><a href="._Bayesian-bs011.html#___sec5" style="font-size: 80%;"><b>Bayes' theorem</b></a></li>
<!-- navigation toc: --> <li><a href="._Bayesian-bs013.html#___sec6" style="font-size: 80%;"><b>The friends of Bayes' theorem</b></a></li>
<!-- navigation toc: --> <li><a href="._Bayesian-bs014.html#___sec7" style="font-size: 80%;"><b>Inference With Parametric Models</b></a></li>
<!-- navigation toc: --> <li><a href="._Bayesian-bs016.html#___sec8" style="font-size: 80%;"><b>Illustrative examples with python code</b></a></li>
<!-- navigation toc: --> <li><a href="._Bayesian-bs017.html#___sec9" style="font-size: 80%;"><b>Example: Is this a fair coin?</b></a></li>
<!-- navigation toc: --> <li><a href="._Bayesian-bs018.html#___sec10" style="font-size: 80%;"><b>A few words on different priors</b></a></li>
<!-- navigation toc: --> <li><a href="._Bayesian-bs019.html#___sec11" style="font-size: 80%;"><b>Bayesian parameter estimation (single parameter)</b></a></li>
<!-- navigation toc: --> <li><a href="._Bayesian-bs020.html#___sec12" style="font-size: 80%;">&nbsp;&nbsp;&nbsp;Example: Measured flux from a star</a></li>
<!-- navigation toc: --> <li><a href="._Bayesian-bs020.html#___sec13" style="font-size: 80%;">&nbsp;&nbsp;&nbsp;Simple Photon Counts: Frequentist Approach</a></li>
<!-- navigation toc: --> <li><a href="._Bayesian-bs020.html#___sec14" style="font-size: 80%;">&nbsp;&nbsp;&nbsp;Simple Photon Counts: Bayesian Approach</a></li>
<!-- navigation toc: --> <li><a href="._Bayesian-bs020.html#___sec15" style="font-size: 80%;">&nbsp;&nbsp;&nbsp;A note about priors</a></li>
<!-- navigation toc: --> <li><a href="._Bayesian-bs020.html#___sec16" style="font-size: 80%;">&nbsp;&nbsp;&nbsp;Simple Photon Counts: Bayesian approach in practice</a></li>
<!-- navigation toc: --> <li><a href="._Bayesian-bs020.html#___sec17" style="font-size: 80%;">&nbsp;&nbsp;&nbsp;Best estimates and confidence intervals</a></li>
<!-- navigation toc: --> <li><a href="._Bayesian-bs020.html#___sec18" style="font-size: 80%;">&nbsp;&nbsp;&nbsp;Simple Photon Counts: Best estimates and confidence intervals</a></li>
<!-- navigation toc: --> <li><a href="._Bayesian-bs021.html#___sec19" style="font-size: 80%;"><b>Bayesian parameter estimation (multiple parameters, covariance)</b></a></li>
<!-- navigation toc: --> <li><a href="._Bayesian-bs022.html#___sec20" style="font-size: 80%;"><b>Bayesian model selection</b></a></li>
</ul>
</li>
@@ -134,7 +174,7 @@ MathJax.Hub.Config({
<center>[3] <b>Department of Physics and Astronomy and National Superconducting Cyclotron Laboratory, Michigan State University</b></center>
<br>
<p>
<center><h4>Mar 10, 2018</h4></center> <!-- date -->
<center><h4>Mar 13, 2018</h4></center> <!-- date -->
<br>
<p>
@@ -158,7 +198,7 @@ MathJax.Hub.Config({
<li><a href="._Bayesian-bs008.html">9</a></li>
<li><a href="._Bayesian-bs009.html">10</a></li>
<li><a href="">...</a></li>
<li><a href="._Bayesian-bs016.html">17</a></li>
<li><a href="._Bayesian-bs022.html">23</a></li>
<li><a href="._Bayesian-bs001.html">&raquo;</a></li>
</ul>
<!-- ------------------- end of main content --------------- -->
doc/pub/Bayesian/html/._Bayesian-bs001.html
+51 -11
View File
@@ -48,7 +48,35 @@ Automatically generated HTML file from DocOnce source
("Bayes' theorem", 2, None, '___sec5'),
("The friends of Bayes' theorem", 2, None, '___sec6'),
('Inference With Parametric Models', 2, None, '___sec7'),
('Illustrative examples with python code', 2, None, '___sec8')]}
('Illustrative examples with python code', 2, None, '___sec8'),
('Example: Is this a fair coin?', 2, None, '___sec9'),
('A few words on different priors', 2, None, '___sec10'),
('Bayesian parameter estimation (single parameter)',
2,
None,
'___sec11'),
('Example: Measured flux from a star', 3, None, '___sec12'),
('Simple Photon Counts: Frequentist Approach',
3,
None,
'___sec13'),
('Simple Photon Counts: Bayesian Approach', 3, None, '___sec14'),
('A note about priors', 3, None, '___sec15'),
('Simple Photon Counts: Bayesian approach in practice',
3,
None,
'___sec16'),
('Best estimates and confidence intervals', 3, None, '___sec17'),
('Simple Photon Counts: Best estimates and confidence intervals',
3,
None,
'___sec18'),
('Bayesian parameter estimation (multiple parameters, '
'covariance)',
2,
None,
'___sec19'),
('Bayesian model selection', 2, None, '___sec20')]}
end of tocinfo -->
<body>
@@ -86,15 +114,27 @@ MathJax.Hub.Config({
<li class="dropdown">
<a href="#" class="dropdown-toggle" data-toggle="dropdown">Contents <b class="caret"></b></a>
<ul class="dropdown-menu">
<!-- navigation toc: --> <li><a href="#___sec0" style="font-size: 80%;">What is Bayesian Statistics</a></li>
<!-- navigation toc: --> <li><a href="._Bayesian-bs002.html#___sec1" style="font-size: 80%;">Inference</a></li>
<!-- navigation toc: --> <li><a href="._Bayesian-bs003.html#___sec2" style="font-size: 80%;">Statistical Inference</a></li>
<!-- navigation toc: --> <li><a href="._Bayesian-bs004.html#___sec3" style="font-size: 80%;">Some history</a></li>
<!-- navigation toc: --> <li><a href="._Bayesian-bs010.html#___sec4" style="font-size: 80%;">The Bayesian recipe</a></li>
<!-- navigation toc: --> <li><a href="._Bayesian-bs011.html#___sec5" style="font-size: 80%;">Bayes' theorem</a></li>
<!-- navigation toc: --> <li><a href="._Bayesian-bs013.html#___sec6" style="font-size: 80%;">The friends of Bayes' theorem</a></li>
<!-- navigation toc: --> <li><a href="._Bayesian-bs014.html#___sec7" style="font-size: 80%;">Inference With Parametric Models</a></li>
<!-- navigation toc: --> <li><a href="._Bayesian-bs016.html#___sec8" style="font-size: 80%;">Illustrative examples with python code</a></li>
<!-- navigation toc: --> <li><a href="#___sec0" style="font-size: 80%;"><b>What is Bayesian Statistics</b></a></li>
<!-- navigation toc: --> <li><a href="._Bayesian-bs002.html#___sec1" style="font-size: 80%;"><b>Inference</b></a></li>
<!-- navigation toc: --> <li><a href="._Bayesian-bs003.html#___sec2" style="font-size: 80%;"><b>Statistical Inference</b></a></li>
<!-- navigation toc: --> <li><a href="._Bayesian-bs004.html#___sec3" style="font-size: 80%;"><b>Some history</b></a></li>
<!-- navigation toc: --> <li><a href="._Bayesian-bs010.html#___sec4" style="font-size: 80%;"><b>The Bayesian recipe</b></a></li>
<!-- navigation toc: --> <li><a href="._Bayesian-bs011.html#___sec5" style="font-size: 80%;"><b>Bayes' theorem</b></a></li>
<!-- navigation toc: --> <li><a href="._Bayesian-bs013.html#___sec6" style="font-size: 80%;"><b>The friends of Bayes' theorem</b></a></li>
<!-- navigation toc: --> <li><a href="._Bayesian-bs014.html#___sec7" style="font-size: 80%;"><b>Inference With Parametric Models</b></a></li>
<!-- navigation toc: --> <li><a href="._Bayesian-bs016.html#___sec8" style="font-size: 80%;"><b>Illustrative examples with python code</b></a></li>
<!-- navigation toc: --> <li><a href="._Bayesian-bs017.html#___sec9" style="font-size: 80%;"><b>Example: Is this a fair coin?</b></a></li>
<!-- navigation toc: --> <li><a href="._Bayesian-bs018.html#___sec10" style="font-size: 80%;"><b>A few words on different priors</b></a></li>
<!-- navigation toc: --> <li><a href="._Bayesian-bs019.html#___sec11" style="font-size: 80%;"><b>Bayesian parameter estimation (single parameter)</b></a></li>
<!-- navigation toc: --> <li><a href="._Bayesian-bs020.html#___sec12" style="font-size: 80%;">&nbsp;&nbsp;&nbsp;Example: Measured flux from a star</a></li>
<!-- navigation toc: --> <li><a href="._Bayesian-bs020.html#___sec13" style="font-size: 80%;">&nbsp;&nbsp;&nbsp;Simple Photon Counts: Frequentist Approach</a></li>
<!-- navigation toc: --> <li><a href="._Bayesian-bs020.html#___sec14" style="font-size: 80%;">&nbsp;&nbsp;&nbsp;Simple Photon Counts: Bayesian Approach</a></li>
<!-- navigation toc: --> <li><a href="._Bayesian-bs020.html#___sec15" style="font-size: 80%;">&nbsp;&nbsp;&nbsp;A note about priors</a></li>
<!-- navigation toc: --> <li><a href="._Bayesian-bs020.html#___sec16" style="font-size: 80%;">&nbsp;&nbsp;&nbsp;Simple Photon Counts: Bayesian approach in practice</a></li>
<!-- navigation toc: --> <li><a href="._Bayesian-bs020.html#___sec17" style="font-size: 80%;">&nbsp;&nbsp;&nbsp;Best estimates and confidence intervals</a></li>
<!-- navigation toc: --> <li><a href="._Bayesian-bs020.html#___sec18" style="font-size: 80%;">&nbsp;&nbsp;&nbsp;Simple Photon Counts: Best estimates and confidence intervals</a></li>
<!-- navigation toc: --> <li><a href="._Bayesian-bs021.html#___sec19" style="font-size: 80%;"><b>Bayesian parameter estimation (multiple parameters, covariance)</b></a></li>
<!-- navigation toc: --> <li><a href="._Bayesian-bs022.html#___sec20" style="font-size: 80%;"><b>Bayesian model selection</b></a></li>
</ul>
</li>
@@ -144,7 +184,7 @@ Morten's original plan: Reminder about probabilities from the statistics section
<li><a href="._Bayesian-bs009.html">10</a></li>
<li><a href="._Bayesian-bs010.html">11</a></li>
<li><a href="">...</a></li>
<li><a href="._Bayesian-bs016.html">17</a></li>
<li><a href="._Bayesian-bs022.html">23</a></li>
<li><a href="._Bayesian-bs002.html">&raquo;</a></li>
</ul>
<!-- ------------------- end of main content --------------- -->
doc/pub/Bayesian/html/Bayesian-bs.html
+52 -12
View File
@@ -48,7 +48,35 @@ Automatically generated HTML file from DocOnce source
("Bayes' theorem", 2, None, '___sec5'),
("The friends of Bayes' theorem", 2, None, '___sec6'),
('Inference With Parametric Models', 2, None, '___sec7'),
('Illustrative examples with python code', 2, None, '___sec8')]}
('Illustrative examples with python code', 2, None, '___sec8'),
('Example: Is this a fair coin?', 2, None, '___sec9'),
('A few words on different priors', 2, None, '___sec10'),
('Bayesian parameter estimation (single parameter)',
2,
None,
'___sec11'),
('Example: Measured flux from a star', 3, None, '___sec12'),
('Simple Photon Counts: Frequentist Approach',
3,
None,
'___sec13'),
('Simple Photon Counts: Bayesian Approach', 3, None, '___sec14'),
('A note about priors', 3, None, '___sec15'),
('Simple Photon Counts: Bayesian approach in practice',
3,
None,
'___sec16'),
('Best estimates and confidence intervals', 3, None, '___sec17'),
('Simple Photon Counts: Best estimates and confidence intervals',
3,
None,
'___sec18'),
('Bayesian parameter estimation (multiple parameters, '
'covariance)',
2,
None,
'___sec19'),
('Bayesian model selection', 2, None, '___sec20')]}
end of tocinfo -->
<body>
@@ -86,15 +114,27 @@ MathJax.Hub.Config({
<li class="dropdown">
<a href="#" class="dropdown-toggle" data-toggle="dropdown">Contents <b class="caret"></b></a>
<ul class="dropdown-menu">
<!-- navigation toc: --> <li><a href="._Bayesian-bs001.html#___sec0" style="font-size: 80%;">What is Bayesian Statistics</a></li>
<!-- navigation toc: --> <li><a href="._Bayesian-bs002.html#___sec1" style="font-size: 80%;">Inference</a></li>
<!-- navigation toc: --> <li><a href="._Bayesian-bs003.html#___sec2" style="font-size: 80%;">Statistical Inference</a></li>
<!-- navigation toc: --> <li><a href="._Bayesian-bs004.html#___sec3" style="font-size: 80%;">Some history</a></li>
<!-- navigation toc: --> <li><a href="._Bayesian-bs010.html#___sec4" style="font-size: 80%;">The Bayesian recipe</a></li>
<!-- navigation toc: --> <li><a href="._Bayesian-bs011.html#___sec5" style="font-size: 80%;">Bayes' theorem</a></li>
<!-- navigation toc: --> <li><a href="._Bayesian-bs013.html#___sec6" style="font-size: 80%;">The friends of Bayes' theorem</a></li>
<!-- navigation toc: --> <li><a href="._Bayesian-bs014.html#___sec7" style="font-size: 80%;">Inference With Parametric Models</a></li>
<!-- navigation toc: --> <li><a href="._Bayesian-bs016.html#___sec8" style="font-size: 80%;">Illustrative examples with python code</a></li>
<!-- navigation toc: --> <li><a href="._Bayesian-bs001.html#___sec0" style="font-size: 80%;"><b>What is Bayesian Statistics</b></a></li>
<!-- navigation toc: --> <li><a href="._Bayesian-bs002.html#___sec1" style="font-size: 80%;"><b>Inference</b></a></li>
<!-- navigation toc: --> <li><a href="._Bayesian-bs003.html#___sec2" style="font-size: 80%;"><b>Statistical Inference</b></a></li>
<!-- navigation toc: --> <li><a href="._Bayesian-bs004.html#___sec3" style="font-size: 80%;"><b>Some history</b></a></li>
<!-- navigation toc: --> <li><a href="._Bayesian-bs010.html#___sec4" style="font-size: 80%;"><b>The Bayesian recipe</b></a></li>
<!-- navigation toc: --> <li><a href="._Bayesian-bs011.html#___sec5" style="font-size: 80%;"><b>Bayes' theorem</b></a></li>
<!-- navigation toc: --> <li><a href="._Bayesian-bs013.html#___sec6" style="font-size: 80%;"><b>The friends of Bayes' theorem</b></a></li>
<!-- navigation toc: --> <li><a href="._Bayesian-bs014.html#___sec7" style="font-size: 80%;"><b>Inference With Parametric Models</b></a></li>
<!-- navigation toc: --> <li><a href="._Bayesian-bs016.html#___sec8" style="font-size: 80%;"><b>Illustrative examples with python code</b></a></li>
<!-- navigation toc: --> <li><a href="._Bayesian-bs017.html#___sec9" style="font-size: 80%;"><b>Example: Is this a fair coin?</b></a></li>
<!-- navigation toc: --> <li><a href="._Bayesian-bs018.html#___sec10" style="font-size: 80%;"><b>A few words on different priors</b></a></li>
<!-- navigation toc: --> <li><a href="._Bayesian-bs019.html#___sec11" style="font-size: 80%;"><b>Bayesian parameter estimation (single parameter)</b></a></li>
<!-- navigation toc: --> <li><a href="._Bayesian-bs020.html#___sec12" style="font-size: 80%;">&nbsp;&nbsp;&nbsp;Example: Measured flux from a star</a></li>
<!-- navigation toc: --> <li><a href="._Bayesian-bs020.html#___sec13" style="font-size: 80%;">&nbsp;&nbsp;&nbsp;Simple Photon Counts: Frequentist Approach</a></li>
<!-- navigation toc: --> <li><a href="._Bayesian-bs020.html#___sec14" style="font-size: 80%;">&nbsp;&nbsp;&nbsp;Simple Photon Counts: Bayesian Approach</a></li>
<!-- navigation toc: --> <li><a href="._Bayesian-bs020.html#___sec15" style="font-size: 80%;">&nbsp;&nbsp;&nbsp;A note about priors</a></li>
<!-- navigation toc: --> <li><a href="._Bayesian-bs020.html#___sec16" style="font-size: 80%;">&nbsp;&nbsp;&nbsp;Simple Photon Counts: Bayesian approach in practice</a></li>
<!-- navigation toc: --> <li><a href="._Bayesian-bs020.html#___sec17" style="font-size: 80%;">&nbsp;&nbsp;&nbsp;Best estimates and confidence intervals</a></li>
<!-- navigation toc: --> <li><a href="._Bayesian-bs020.html#___sec18" style="font-size: 80%;">&nbsp;&nbsp;&nbsp;Simple Photon Counts: Best estimates and confidence intervals</a></li>
<!-- navigation toc: --> <li><a href="._Bayesian-bs021.html#___sec19" style="font-size: 80%;"><b>Bayesian parameter estimation (multiple parameters, covariance)</b></a></li>
<!-- navigation toc: --> <li><a href="._Bayesian-bs022.html#___sec20" style="font-size: 80%;"><b>Bayesian model selection</b></a></li>
</ul>
</li>
@@ -134,7 +174,7 @@ MathJax.Hub.Config({
<center>[3] <b>Department of Physics and Astronomy and National Superconducting Cyclotron Laboratory, Michigan State University</b></center>
<br>
<p>
<center><h4>Mar 10, 2018</h4></center> <!-- date -->
<center><h4>Mar 13, 2018</h4></center> <!-- date -->
<br>
<p>
@@ -158,7 +198,7 @@ MathJax.Hub.Config({
<li><a href="._Bayesian-bs008.html">9</a></li>
<li><a href="._Bayesian-bs009.html">10</a></li>
<li><a href="">...</a></li>
<li><a href="._Bayesian-bs016.html">17</a></li>
<li><a href="._Bayesian-bs022.html">23</a></li>
<li><a href="._Bayesian-bs001.html">&raquo;</a></li>
</ul>
<!-- ------------------- end of main content --------------- -->
doc/pub/Bayesian/html/Bayesian-reveal.html
+486 -6
View File
@@ -153,7 +153,7 @@ MathJax.Hub.Config({
<center>[3] <b>Department of Physics and Astronomy and National Superconducting Cyclotron Laboratory, Michigan State University</b></center>
<br>
<p>&nbsp;<br>
<center><h4>Mar 10, 2018</h4></center> <!-- date -->
<center><h4>Mar 13, 2018</h4></center> <!-- date -->
<br>
<p>
@@ -314,11 +314,10 @@ $$
The importance of this property to data analysis becomes apparent if we replace \( A \) and \( B \) by hypothesis(\( H \)) and data(\( D \)):
<p>&nbsp;<br>
$$
$$
<p>&nbsp;<br>
p(H|D,I) = \frac{p(D|H,I) p(H|I)}{p(D|I)}.
<p>&nbsp;<br>
$$
\begin{align}
p(H|D,I) &= \frac{p(D|H,I) p(H|I)}{p(D|I)}.
\tag{1}
\end{align}
$$
<p>&nbsp;<br>
@@ -453,6 +452,487 @@ $$
</section>
<section>
<h2 id="___sec9">Example: Is this a fair coin? </h2>
Let us begin with the analysis of data from a simple coin-tossing experiment.
Given that we had observed 6 heads in 8 flips, would you think it was a fair coin? By fair, we mean that we would be prepared to lay an even 1 : 1 bet on the outcome of a flip being a head or a tail. If we decide that the coin was fair, the question which follows naturally is how sure are we that this was so; if it was not fair, how unfair do we think it was? Furthermore, if we were to continue collecting data for this particular coin, observing the outcomes of additional flips, how would we update our belief on the fairness of the coin?
<p>
A sensible way of formulating this problem is to consider a large number of hypotheses about the range in which the bias-weighting of the coin might lie. If we denote the bias-weighting by \( H \), then \( H = 0 \) and \( H = 1 \) can represent a coin which produces a tail or a head on every flip, respectively. There is a continuum of possibilities for the value of H between these limits, with \( H = 0.5 \) indicating a fair coin. Our state of knowledge about the fairness, or the degree of unfairness, of the coin is then completely summarized by specifying how much we believe these various propositions to be true.
<p>
Let us perform a computer simulation of a coin-tossing experiment. This provides the data that we will be analysing.
<p>
<!-- code=python (!bc pycod) typeset with pygments style "perldoc" -->
<div class="highlight" style="background: #eeeedd"><pre style="font-size: 80%; line-height: 125%"><span></span><span style="color: #8B008B; font-weight: bold">import</span> <span style="color: #008b45; text-decoration: underline">numpy</span> <span style="color: #8B008B; font-weight: bold">as</span> <span style="color: #008b45; text-decoration: underline">np</span>
<span style="color: #8B008B; font-weight: bold">import</span> <span style="color: #008b45; text-decoration: underline">matplotlib.pyplot</span> <span style="color: #8B008B; font-weight: bold">as</span> <span style="color: #008b45; text-decoration: underline">plt</span>
</pre></div>
<p>
<!-- code=python (!bc pycod) typeset with pygments style "perldoc" -->
<div class="highlight" style="background: #eeeedd"><pre style="font-size: 80%; line-height: 125%"><span></span>np.random.seed(<span style="color: #B452CD">999</span>) <span style="color: #228B22"># for reproducibility</span>
a=<span style="color: #B452CD">0.6</span> <span style="color: #228B22"># biased coin</span>
flips=np.random.rand(<span style="color: #B452CD">2</span>**<span style="color: #B452CD">12</span>) <span style="color: #228B22"># simulates 4096 coin flips</span>
heads=flips&lt;a <span style="color: #228B22"># boolean array, heads[i]=True if flip i is heads</span>
</pre></div>
<p>
In the light of this data, our inference about the fairness of this coin is summarized by the conditional pdf: \( p(H|D,I) \). This is, of course, shorthand for the limiting case of a continuum of propositions for the value of \( H \); that is to say, the probability that \( H \) lies in an infinitesimally narrow range is given by \( p(H|D,I) dH \).
<p>
To estimate this posterior pdf, we need to use Bayes&#8217; theorem <a href="#mjx-eqn-1">(1)</a>. We will ignore the denominator \( p(D|I) \) as it does not involve bias-weighting explicitly, and it will therefore not affect the shape of the desired pdf. At the end we can evaluate the missing constant subsequently from the normalization condition
<p>&nbsp;<br>
$$
\begin{equation}
\int_0^1 p(H|D,I) dH = 1.
\tag{2}
\end{equation}
$$
<p>&nbsp;<br>
<p>
The prior pdf, \( p(H|I) \), represents what we know about the coin given only the information \( I \) that we are dealing with a &#8216;strange coin&#8217;. We could keep a very open mind about the nature of the coin; a simple probability assignment which reflects this is a uniform, or flat, prior
<p>&nbsp;<br>
$$
\begin{equation}
p(H|I) = \left\{ \begin{array}{ll}
1 & 0 \le H \le 1, \\
0 & \mathrm{otherwise}.
\end{array} \right.
\tag{3}
\end{equation}
$$
<p>&nbsp;<br>
We will get back later to the choice of prior and its effect on the analysis.
<p>
This prior state of knowledge, or ignorance, is modified by the data through the likelihood function \( p(D|H,I) \). It is a measure of the chance that we would have obtained the data that we actually observed, if the value of the bias-weighting was given (as known). If, in the conditioning information \( I \), we assume that the flips of the coin were independent events, so that the outcome of one did not influence that of another, then the probability of obtaining the data `R heads in N tosses' is given by the binomial distribution (we leave a formal definition of this to a statistics textbook)
<p>&nbsp;<br>
$$
\begin{equation}
p(D|H,I) \propto H^R (1-H)^{N-R}.
\tag{4}
\end{equation}
$$
<p>&nbsp;<br>
It seems reasonable because \( H \) is the chance of obtaining a head on any flip, and there were \( R \) of them, and \( 1-H \) is the corresponding probability for a tail, of which there were \( N-R \). We note that this binomial distribution also contains a normalization factor, but we will ignore it since it does not depend explicitly on \( H \), the quantity of interest. It will be absorbed by the normalization condition <a href="#mjx-eqn-2">(2)</a>.
<p>
We perform the setup of this Bayesian framework on the computer.
<p>
<!-- code=python (!bc pycod) typeset with pygments style "perldoc" -->
<div class="highlight" style="background: #eeeedd"><pre style="font-size: 80%; line-height: 125%"><span></span><span style="color: #8B008B; font-weight: bold">def</span> <span style="color: #008b45">prior</span>(H):
p=np.zeros_like(H)
p[(<span style="color: #B452CD">0</span>&lt;=x)&amp;(x&lt;=<span style="color: #B452CD">1</span>)]=<span style="color: #B452CD">1</span> <span style="color: #228B22"># allowed range: 0&lt;=H&lt;=1</span>
<span style="color: #8B008B; font-weight: bold">return</span> p <span style="color: #228B22"># uniform prior</span>
<span style="color: #8B008B; font-weight: bold">def</span> <span style="color: #008b45">likelihood</span>(H,data):
N = <span style="color: #658b00">len</span>(data)
no_of_heads = <span style="color: #658b00">sum</span>(data)
no_of_tails = N - no_of_heads
<span style="color: #8B008B; font-weight: bold">return</span> H**no_of_heads * (<span style="color: #B452CD">1</span>-H)**no_of_tails
<span style="color: #8B008B; font-weight: bold">def</span> <span style="color: #008b45">posterior</span>(H,data):
p=prior(H)*likelihood(H,data)
norm=np.trapz(p,H)
<span style="color: #8B008B; font-weight: bold">return</span> p/norm
</pre></div>
<p>
The next step is to confront this setup with the simulated data. To get a feel for the result, it is instructive to see how the posterior pdf evolves as we obtain more and more data pertaining to the coin. The results of such an analyses is shown in Fig. <a href="#fig:coinflipping">1</a>.
<p>
<!-- code=python (!bc pycod) typeset with pygments style "perldoc" -->
<div class="highlight" style="background: #eeeedd"><pre style="font-size: 80%; line-height: 125%"><span></span>x=np.linspace(<span style="color: #B452CD">0</span>,<span style="color: #B452CD">1</span>,<span style="color: #B452CD">100</span>)
fig, axs = plt.subplots(nrows=<span style="color: #B452CD">4</span>,ncols=<span style="color: #B452CD">3</span>,sharex=<span style="color: #658b00">True</span>,sharey=<span style="color: #CD5555">&#39;row&#39;</span>)
axs_vec=np.reshape(axs,-<span style="color: #B452CD">1</span>)
axs_vec[<span style="color: #B452CD">0</span>].plot(x,prior(x))
<span style="color: #8B008B; font-weight: bold">for</span> ndouble <span style="color: #8B008B">in</span> <span style="color: #658b00">range</span>(<span style="color: #B452CD">11</span>):
ax=axs_vec[<span style="color: #B452CD">1</span>+ndouble]
ax.plot(x,posterior(x,heads[:<span style="color: #B452CD">2</span>**ndouble]))
ax.text(<span style="color: #B452CD">0.1</span>, <span style="color: #B452CD">0.8</span>, <span style="color: #CD5555">&#39;$N={0}$&#39;</span>.format(<span style="color: #B452CD">2</span>**ndouble), transform=ax.transAxes)
<span style="color: #8B008B; font-weight: bold">for</span> row <span style="color: #8B008B">in</span> <span style="color: #658b00">range</span>(<span style="color: #B452CD">4</span>): axs[row,<span style="color: #B452CD">0</span>].set_ylabel(<span style="color: #CD5555">&#39;$p(H|D_\mathrm{obs},I)$&#39;</span>)
<span style="color: #8B008B; font-weight: bold">for</span> col <span style="color: #8B008B">in</span> <span style="color: #658b00">range</span>(<span style="color: #B452CD">3</span>): axs[-<span style="color: #B452CD">1</span>,col].set_xlabel(<span style="color: #CD5555">&#39;$H$&#39;</span>)
</pre></div>
<p>
<center> <!-- figure label: --> <div id="fig:coinflipping"></div> <!-- FIGURE -->
<hr class="figure">
<center><p class="caption">Figure 1: The evolution of the posterior pdf for the bias-weighting of a coin, as the number of data available increases. The figure on the top left-hand corner of each panel shows the number of data included in the analysis. <!-- caption label: fig:coinflipping --> </p></center>
<p><img src="fig/coinflipping_fig_1.png" align="bottom" width=500></p>
</center>
<p>
The panel in the top left-hand corner shows the posterior pdf for \( H \) given no data, i.e., it is the same as the prior pdf of Eq. <a href="#mjx-eqn-3">(3)</a>. It indicates that we have no more reason to believe that the coin is fair than we have to think that it is double-headed, double-tailed, or of any other intermediate bias-weighting.
<p>
The first flip is obviously tails. At this point we have no evidence that the coin has a side with heads, as indicated by the pdf going to zero as \( H \to 1 \). The second flip is obviously heads and we have now excluded both extreme options \( H=0 \) (double-tailed) and \( H=1 \) (double-headed). We can note that the posterior at this point has the simple form \( p(H|D,I) = H(1-H) \) for \( 0 \le H \le 1 \).
<p>
The remainder of Fig. <a href="#fig:coinflipping">1</a> shows how the posterior pdf evolves as the number of data analysed becomes larger and larger. We see that the position of the maximum moves around, but that the amount by which it does so decreases with the increasing number of observations. The width of the posterior pdf also becomes narrower with more data, indicating that we are becoming increasingly confident in our estimate of the bias-weighting. For the coin in this example, the best estimate of \( H \) eventually converges to 0.6, which, of course, was the value chosen to simulate the flips.
</section>
<section>
<h2 id="___sec10">A few words on different priors </h2>
<ul>
<p><li> uniform</li>
<p><li> Gaussian</li>
<p><li> Jeffrey's prior</li>
</ul>
<p>
Repeat the coin flipping experiment with other priors.
</section>
<section>
<h2 id="___sec11">Bayesian parameter estimation (single parameter) </h2>
<div class="alert alert-block alert-block alert-text-normal">
<b></b>
<p>
We will now consider the very important task of model parameter estimation using statistical inference.
<!-- begin inline comment -->
<font color="red">(<b>CF 1</b>: maybe stress that model parameters are not random variables, and the meaning of parameter estimation is therefore very different between frequentist and bayesian approaches.)</font>
<!-- end inline comment -->
<p>
Throughout this section we will consider a specific example that involves a model with a single parameter: &quot;Measured flux from a star&quot;.
</div>
</section>
<section>
<h3 id="___sec12">Example: Measured flux from a star </h3>
Adapted from the blog <a href="http://jakevdp.github.io" target="_blank">Pythonic Perambulations</a> by Jake VanderPlas.
<p>
Imagine that we point our telescope to the sky, and observe the light coming from a single star. For the time being, we'll assume that the star's true flux is constant with time, i.e. that is it has a fixed value \( F_\mathrm{true} \) (we'll also ignore effects like sky noise and other sources of systematic error). We'll assume that we perform a series of \( N \) measurements with our telescope, where the ith measurement reports the observed photon flux \( F_i \) and error \( e_i \) [<a id="link_footnote_2" href="#def_footnote_2">2</a>].
The question is, given this set of measurements \( D = \{F_i, e_i\} \), what is our best estimate of the true flux \( F_\mathrm{true} \)?
<p id="def_footnote_2"><a href="#link_footnote_2"><b>2:</b></a> We'll make the reasonable assumption that errors are Gaussian. In a Frequentist perspective, \( e_i \) is the standard deviation of the results of a single measurement event in the limit of repetitions of <em>that event</em>. In the Bayesian perspective, \( e_i \) is the standard deviation of the (Gaussian) probability distribution describing our knowledge of that particular measurement given its observed value.</p>
<p>
Because the measurements are number counts, a Poisson distribution is a good approximation to the measurement process:
<p>
<!-- code=python (!bc pycod) typeset with pygments style "perldoc" -->
<div class="highlight" style="background: #eeeedd"><pre style="font-size: 80%; line-height: 125%"><span></span>np.random.seed(<span style="color: #B452CD">1</span>) <span style="color: #228B22"># for repeatability</span>
F_true = <span style="color: #B452CD">1000</span> <span style="color: #228B22"># true flux, say number of photons measured in 1 second</span>
N = <span style="color: #B452CD">50</span> <span style="color: #228B22"># number of measurements</span>
F = stats.poisson(F_true).rvs(N)
<span style="color: #228B22"># N measurements of the flux</span>
e = np.sqrt(F) <span style="color: #228B22"># errors on Poisson counts estimated via square root</span>
</pre></div>
<p>
Now let's make a simple visualization of the &quot;observed&quot; data, see Fig. <a href="#fig:flux">2</a>.
<p>
<!-- code=python (!bc pycod) typeset with pygments style "perldoc" -->
<div class="highlight" style="background: #eeeedd"><pre style="font-size: 80%; line-height: 125%"><span></span>fig, ax = plt.subplots()
ax.errorbar(F, np.arange(N), xerr=e, fmt=<span style="color: #CD5555">&#39;ok&#39;</span>, ecolor=<span style="color: #CD5555">&#39;gray&#39;</span>, alpha=<span style="color: #B452CD">0.5</span>)
ax.vlines([F_true], <span style="color: #B452CD">0</span>, N, linewidth=<span style="color: #B452CD">5</span>, alpha=<span style="color: #B452CD">0.2</span>)
ax.set_xlabel(<span style="color: #CD5555">&quot;Flux&quot;</span>);ax.set_ylabel(<span style="color: #CD5555">&quot;measurement number&quot;</span>);
</pre></div>
<p>
<center> <!-- figure label: --> <div id="fig:flux"></div> <!-- FIGURE -->
<hr class="figure">
<center><p class="caption">Figure 2: Single photon counts (flux measurements). <!-- caption label: fig:flux --> </p></center>
<p><img src="fig/singlephotoncount_fig_1.png" align="bottom" width=400></p>
</center>
<p>
These measurements each have a different error \( e_i \) which is estimated from Poisson statistics using the standard square-root rule. In this toy example we already know the true flux \( F_\mathrm{true} \), but the question is this: given our measurements and errors, what is our best estimate of the true flux?
<p>
Let's take a look at the frequentist and Bayesian approaches to solving this.
<h3 id="___sec13">Simple Photon Counts: Frequentist Approach </h3>
We'll start with the classical frequentist maximum likelihood approach. Given a single observation \( D_i = (F_i, e_i) \), we can compute the probability distribution of the measurement given the true flux Ftrue given our assumption of Gaussian errors
<p>&nbsp;<br>
$$
\begin{equation}
p(D_i | F_\mathrm{true}, I) = \frac{1}{\sqrt{2\pi e_i^2}} \exp \left( \frac{-(F_i-F_\mathrm{true})^2}{2e_i^2} \right).
\tag{5}
\end{equation}
$$
<p>&nbsp;<br>
This should be read &quot;the probability of \( D_i \) given \( F_\mathrm{true} \)
equals ...&quot;. You should recognize this as a normal distribution with mean \( F_\mathrm{true} \) and standard deviation \( e_i \).
<p>
We construct the <em>likelihood function</em> by computing the product of the probabilities for each data point
<p>&nbsp;<br>
$$
\begin{equation}
\mathcal{L}(D | F_\mathrm{true}, I) = \prod_{i=1}^N p(D_i | F_\mathrm{true}, I),
\tag{6}
\end{equation}
$$
<p>&nbsp;<br>
here \( D = \{D_i\} \) represents the entire set of measurements. Because the value of the likelihood can become very small, it is often more convenient to instead compute the log-likelihood. Combining the previous two equations and computing the log, we have
<p>&nbsp;<br>
$$
\begin{equation}
\log\mathcal{L} = -\frac{1}{2} \sum_{i=1}^N \left[ \log(2\pi e_i^2) + \frac{(F_i-F_\mathrm{true})^2}{e_i^2} \right].
\tag{7}
\end{equation}
$$
<p>&nbsp;<br>
<p>
What we'd like to do is determine \( F_\mathrm{true} \) such that the likelihood is maximized. For this simple problem, the maximization can be computed analytically (i.e. by setting \( d\log\mathcal{L}/d F_\mathrm{true} = 0 \)). This results in the following observed estimate of \( F_\mathrm{true} \)
<p>&nbsp;<br>
$$
\begin{equation}
F_\mathrm{est} = \sum_{i=1}^N w_i F_i; \quad w_i = 1/e_i^2.
\tag{8}
\end{equation}
$$
<p>&nbsp;<br>
Notice that in the special case of all errors \( e_i \) being equal, this reduces to
<p>&nbsp;<br>
$$
\begin{equation}
F_\mathrm{est} = \frac{1}{N} \sum_{i=1} F_i.
\tag{9}
\end{equation}
$$
<p>&nbsp;<br>
That is, in agreement with intuition, \( F_\mathrm{est} \) is simply the mean of the observed data when errors are equal.
<p>
We can go further and ask what the error of our estimate is. In the frequentist approach, this can be accomplished by fitting a Gaussian approximation to the likelihood curve at maximum; in this simple case this can also be solved analytically (the sum of Gaussians is also a Gaussian). It can be shown that the standard deviation of this Gaussian approximation is
<p>&nbsp;<br>
$$
\begin{equation}
\sigma_\mathrm{est} = \sum_{i=1}^N w_i.
\tag{10}
\end{equation}
$$
<p>&nbsp;<br>
These results are fairly simple calculations; let's evaluate them for our toy dataset:
<p>
<!-- code=python (!bc pycod) typeset with pygments style "perldoc" -->
<div class="highlight" style="background: #eeeedd"><pre style="font-size: 80%; line-height: 125%"><span></span>w=<span style="color: #B452CD">1.</span>/e**<span style="color: #B452CD">2</span>
<span style="color: #8B008B; font-weight: bold">print</span>(<span style="color: #CD5555">&quot;&quot;&quot;</span>
<span style="color: #CD5555">F_true = {0}</span>
<span style="color: #CD5555">F_est = {1:.0f} +/- {2:.0f} (based on {3} measurements) &quot;&quot;&quot;</span>\
.format(F_true, (w * F).sum() / w.sum(), w.sum() ** -<span style="color: #B452CD">0.5</span>, N))
</pre></div>
<p>
<code>F_true = 1000</code> <br />
<code>F_est = 998 +/- 4 (based on 50 measurements)</code> <br />
<p>
We find that for 50 measurements of the flux, our estimate has an error of about 0.4% and is consistent with the input value.
<h3 id="___sec14">Simple Photon Counts: Bayesian Approach </h3>
The Bayesian approach, as you might expect, begins and ends with probabilities. Our hypothesis is that the star has a constant flux \( F_\mathrm{true} \). It recognizes that what we fundamentally want to compute is our knowledge of the parameters in question given the data and other information (such as our knowledge of uncertainties for the observed values), i.e. in this case, \( p(F_\mathrm{true} | D,I) \).
Note that this formulation of the problem is fundamentally contrary to the frequentist philosophy, which says that probabilities have no meaning for model parameters like \( F_\mathrm{true} \). Nevertheless, within the Bayesian philosophy this is perfectly acceptable.
<p>
To compute this result, Bayesians next apply Bayes' Theorem <a href="#mjx-eqn-1">(1)</a>.
If we set the prior \( p(F_\mathrm{true}|I) \propto 1 \) (a flat prior), we find
\( p(F_\mathrm{true}|D,I) \propto p(D | F_\mathrm{true},I) \equiv \mathcal{L}(D | F_\mathrm{true},I) \)
and the Bayesian probability is maximized at precisely the same value as the frequentist result! So despite the philosophical differences, we see that (for this simple problem at least) the Bayesian and frequentist point estimates are equivalent.
<h3 id="___sec15">A note about priors </h3>
The prior allows inclusion of other information into the computation, which becomes very useful in cases where multiple measurement strategies are being combined to constrain a single model. The necessity to specify a prior, however, is one of the more controversial pieces of Bayesian analysis.
A frequentist will point out that the prior is problematic when no true prior information is available. Though it might seem straightforward to use a noninformative prior like the flat prior mentioned above, there are some <a href="http://normaldeviate.wordpress.com/2013/07/13/lost-causes-in-statistics-ii-noninformative- priors/comment-page-1/" target="_blank">surprisingly subtleties</a> involved. It turns out that in many situations, a truly noninformative prior does not exist! Frequentists point out that the subjective choice of a prior which necessarily biases your result has no place in statistical data analysis.
A Bayesian would counter that frequentism doesn't solve this problem, but simply skirts the question. Frequentism can often be viewed as simply a special case of the Bayesian approach for some (implicit) choice of the prior: a Bayesian would say that it's better to make this implicit choice explicit, even if the choice might include some subjectivity.
<h3 id="___sec16">Simple Photon Counts: Bayesian approach in practice </h3>
Leaving these philosophical debates aside for the time being, let's address how Bayesian results are generally computed in practice. For a one parameter problem like the one considered here, it's as simple as computing the posterior probability \( p(F_\mathrm{true} | D,I) \) as a function of \( F_\mathrm{true} \): this is the distribution reflecting our knowledge of the parameter \( F_\mathrm{true} \).
But as the dimension of the model grows, this direct approach becomes increasingly intractable. For this reason, Bayesian calculations often depend on sampling methods such as Markov Chain Monte Carlo (MCMC). For this practical example, let us apply an MCMC approach using Dan Foreman-Mackey's <a href="http://dan.iel.fm/emcee/current/" target="_blank">emcee</a> package. Keep in mind here that the goal is to generate a set of points drawn from the posterior probability distribution, and to use those points to determine the answer we seek.
To perform this MCMC, we start by defining Python functions for the prior \( p(F_\mathrm{true} | I) \), the likelihood \( p(D | F_\mathrm{true},I) \), and the posterior \( p(F_\mathrm{true} | D,I) \), noting that none of these need be properly normalized. Our model here is one-dimensional, but to handle multi-dimensional models we'll define the model in terms of an array of parameters \( \boldsymbol{\alpha} \), which in this case is \( \boldsymbol{\alpha} = [F_\mathrm{true}] \)
<p>
<!-- code=python (!bc pycod) typeset with pygments style "perldoc" -->
<div class="highlight" style="background: #eeeedd"><pre style="font-size: 80%; line-height: 125%"><span></span><span style="color: #8B008B; font-weight: bold">def</span> <span style="color: #008b45">log_prior</span>(alpha):
<span style="color: #8B008B; font-weight: bold">return</span> <span style="color: #B452CD">0</span> <span style="color: #228B22"># flat prior</span>
<span style="color: #8B008B; font-weight: bold">def</span> <span style="color: #008b45">log_likelihood</span>(alpha, F, e):
<span style="color: #8B008B; font-weight: bold">return</span> -<span style="color: #B452CD">0.5</span> * np.sum(np.log(<span style="color: #B452CD">2</span> * np.pi * e ** <span style="color: #B452CD">2</span>) \
+ (F - alpha[<span style="color: #B452CD">0</span>]) ** <span style="color: #B452CD">2</span> / e ** <span style="color: #B452CD">2</span>)
<span style="color: #8B008B; font-weight: bold">def</span> <span style="color: #008b45">log_posterior</span>(alpha, F, e):
<span style="color: #8B008B; font-weight: bold">return</span> log_prior(alpha) + log_likelihood(alpha, F, e)
</pre></div>
<p>
Now we set up the problem, including generating some random starting guesses for the multiple chains of points.
<p>
<!-- code=python (!bc pycod) typeset with pygments style "perldoc" -->
<div class="highlight" style="background: #eeeedd"><pre style="font-size: 80%; line-height: 125%"><span></span>ndim = <span style="color: #B452CD">1</span> <span style="color: #228B22"># number of parameters in the model</span>
nwalkers = <span style="color: #B452CD">50</span> <span style="color: #228B22"># number of MCMC walkers</span>
nburn = <span style="color: #B452CD">1000</span> <span style="color: #228B22"># &quot;burn-in&quot; period to let chains stabilize</span>
nsteps = <span style="color: #B452CD">2000</span> <span style="color: #228B22"># number of MCMC steps to take</span>
<span style="color: #228B22"># we&#39;ll start at random locations between 0 and 2000</span>
starting_guesses = <span style="color: #B452CD">2000</span> * np.random.rand(nwalkers, ndim)
sampler = emcee.EnsembleSampler(nwalkers, ndim, log_posterior, args=[F,e])
sampler.run_mcmc(starting_guesses, nsteps)
<span style="color: #228B22"># Shape of sampler.chain = (nwalkers, nsteps, ndim)</span>
<span style="color: #228B22"># Flatten the sampler chain and discard burn-in points:</span>
samples = sampler.chain[:, nburn:, :].reshape((-<span style="color: #B452CD">1</span>, ndim))
</pre></div>
<p>
If this all worked correctly, the array sample should contain a series of 50,000 points drawn from the posterior. Let's plot them and check. See results in Fig. <a href="#fig:flux-bayesian">3</a>.
<p>
<!-- code=python (!bc pycod) typeset with pygments style "perldoc" -->
<div class="highlight" style="background: #eeeedd"><pre style="font-size: 80%; line-height: 125%"><span></span>fig, ax = plt.subplots()
ax.hist(samples, bins=<span style="color: #B452CD">50</span>, histtype=<span style="color: #CD5555">&quot;stepfilled&quot;</span>, alpha=<span style="color: #B452CD">0.3</span>, normed=<span style="color: #658b00">True</span>)
ax.set_xlabel(<span style="color: #CD5555">r&#39;$F_\mathrm{est}$&#39;</span>)
ax.set_ylabel(<span style="color: #CD5555">r&#39;$p(F_\mathrm{est}|D,I)$&#39;</span>)
</pre></div>
<p>
<center> <!-- figure label: --> <div id="fig:flux-bayesian"></div> <!-- FIGURE -->
<hr class="figure">
<center><p class="caption">Figure 3: Bayesian posterior pdf (represented by a histogram of MCMC samples) from flux measurements. <!-- caption label: fig:flux-bayesian --> </p></center>
<p><img src="fig/singlephotoncount_fig_2.png" align="bottom" width=400></p>
</center>
<h3 id="___sec17">Best estimates and confidence intervals </h3>
The posterior distribution from our Bayesian data analysis is the key quantity that encodes our inference about the values of the model parameters, given the data and the relevant background information. Often, however, we wish to summarize this result with just a few numbers: the best estimate and a measure of its reliability.
<p>
There are a few different options for this. The choice of the most appropriate one depends mainly on the shape of the posterior distribution:
<p>
<em>Symmetric posterior pdfs</em>: Since the probability (density) associated with any particular value of the parameter is a measure of how much we believe that it lies in the neighbourhood of that point, our best estimate is given by the maximum of the posterior pdf. If we denote the quantity of interest by \( X \), with a posterior pdf \( P =p(X|D,I) \), then the best estimate of its value \( X_0 \) is given by the condition \( dP/dX|_{X=X_0}=0 \). Strictly speaking, we should also check the sign of the second derivative to ensure that \( X_0 \) represents a maximum.
<p>
To obtain a measure of the reliability of this best estimate, we need to look at the width or spread of the posterior pdf about \( X_0 \). When considering the behaviour of any function in the neighbourhood of a particular point, it is often helpful to carry out a Taylor series expansion; this is simply a standard tool for (locally) approximating a complicated function by a low-order polynomial. The linear term is zero at the maximum and the quadratic term is often the dominating one determining the width of the posterior pdf. Ignoring all the higher-order terms we arrive at the Gaussian approximation
<p>&nbsp;<br>
$$
\begin{equation}
p(X|D,I) \approx \frac{1}{\sigma\sqrt{2\pi}} \exp \left[ -\frac{(x-\mu)^2}{2\sigma^2} \right],
\tag{11}
\end{equation}
$$
<p>&nbsp;<br>
where the mean \( \mu = X_0 \) and the variance \( \sigma = \left( - \left. \frac{d^2L}{dX^2} \right|_{X_0} \right)^{-1/2} \), where \( L \) is the logarithm of the posterior \( P \). Our inference about the quantity of interest is conveyed very concisely, therefore, by the statement \( X = X_0 \pm \sigma \), and
<p>&nbsp;<br>
$$
$$
<p>&nbsp;<br>
p(X_0-\sigma < X < X_0+\sigma | D,I) = \int_{X_0-\sigma}^{X_0+\sigma} p(X|D,I) dX \approx 0.67.
<p>&nbsp;<br>
$$
$$
<p>&nbsp;<br>
<p>
<em>Asymmetric posterior pdfs</em>: While the maximum of the posterior (\( X_0 \)) can still be regarded as giving the best estimate, the true value is now more likely to be on one side of this rather than the other. Alternatively one can compute the mean value, \( \langle X \rangle = \int X p(X|D,I) dX \), although this tends to overemphasise very long tails. The best option is probably a compromise that can be employed when having access to a large sample from the posterior (as provided by an MCMC), namely to give the median of this ensamble.
<p>
Furthermore, the concept of an error-bar does not seem appropriate in this case, as it implicitly entails the idea of symmetry. A good way of expressing the reliability with which a parameter can be inferred, for an asymmetric posterior pdf, is rather through a <em>confidence interval</em>. Since the area under the posterior pdf between \( X_1 \) and \( X_2 \) is proportional to how much we believe that \( X \) lies in that range, the shortest interval that encloses 67% of the area represents a sensible measure of the uncertainty of the estimate. Obviously we can choose to provide some other degree-of-belief that we think is relevant for the case at hand. Assuming that the posterior pdf has been normalized, to have unit area, we need to find \( X_1 \) and \( X_2 \) such that:
<p>&nbsp;<br>
$$
$$
<p>&nbsp;<br>
p(X_1 < X < X_2 | D,I) = \int_{X_1}^{X_2} p(X|D,I) dX \approx 0.67,
<p>&nbsp;<br>
$$
$$
<p>&nbsp;<br>
where the difference \( X_2 - X_1 \) is as small as possible. The region \( X_1 < X < X_2 \) is then called the shortest 67% confidence interval.
<p>
<em>Multimodal posterior pdfs</em>: We can sometimes obtain posteriors which are multimodal; i.e. contains several disconnected regions with large probabilities. There is no difficulty when one of the maxima is very much larger than the others: we can simply ignore the subsidiary solutions, to a good approximation, and concentrate on the global maximum. The problem arises when there are several maxima of comparable magnitude. What do we now mean by a best estimate, and how should we quantify its reliability? The idea of a best estimate and an error-bar, or even a confidence interval, is merely an attempt to summarize the posterior with just two or three numbers; sometimes this just can&#8217;t be done, and so these concepts are not valid. For the bimodal case we might be able to characterize the posterior in terms of a few numbers: two best estimates and their associated error-bars, or disjoint confidence intervals. For a general multimodal pdf, the most honest thing we can do is just display the posterior itself.
<h3 id="___sec18">Simple Photon Counts: Best estimates and confidence intervals </h3>
To compute these numbers for our example, you would run:
<p>
<!-- code=python (!bc pycod) typeset with pygments style "perldoc" -->
<div class="highlight" style="background: #eeeedd"><pre style="font-size: 80%; line-height: 125%"><span></span>sampper=np.percentile(samples, [<span style="color: #B452CD">2.5</span>, <span style="color: #B452CD">16.5</span>, <span style="color: #B452CD">50</span>, <span style="color: #B452CD">83.5</span>, <span style="color: #B452CD">97.5</span>],axis=<span style="color: #B452CD">0</span>).flatten()
<span style="color: #8B008B; font-weight: bold">print</span>(<span style="color: #CD5555">&quot;&quot;&quot;</span>
<span style="color: #CD5555">F_true = {0}</span>
<span style="color: #CD5555">Based on {1} measurements the posterior point estimates are:</span>
<span style="color: #CD5555">...F_est = {2:.0f} +/- {3:.0f}</span>
<span style="color: #CD5555">or using credible intervals:</span>
<span style="color: #CD5555">...F_est = {4:.0f} (posterior median) </span>
<span style="color: #CD5555">...F_est in [{5:.0f}, {6:.0f}] (67% credible interval) </span>
<span style="color: #CD5555">...F_est in [{7:.0f}, {8:.0f}] (95% credible interval) &quot;&quot;&quot;</span>\
.format(F_true, N, np.mean(samples), np.std(samples), \
sampper[<span style="color: #B452CD">2</span>], sampper[<span style="color: #B452CD">1</span>], sampper[<span style="color: #B452CD">3</span>], sampper[<span style="color: #B452CD">0</span>], sampper[<span style="color: #B452CD">4</span>]))
</pre></div>
<p>
<code>F_true = 1000</code> <br />
<code>Based on 50 measurements the posterior point estimates are:</code> <br />
<code>...F_est = 998 +/- 4</code> <br />
<code>or using credible intervals:</code> <br />
<code>...F_est = 998 (posterior median)</code> <br />
<code>...F_est in [993, 1002] (67% credible interval)</code> <br />
<code>...F_est in [989, 1006] (95% credible interval)</code> <br />
<p>
In this particular example, the posterior pdf is actually a Gaussian (since it is constructed as a product of Gaussians), and the mean and variance from the quadratic approximation will agree exactly with the frequentist approach.
<p>
From this final result you might come away with the impression that the Bayesian method is unnecessarily complicated, and in this case it certainly is. Using an MCMC sampler to characterize a one-dimensional normal distribution is a bit like using the Death Star to destroy a beach ball, but we did this here because it demonstrates an approach that can scale to complicated posteriors in many, many dimensions, and can provide nice results in more complicated situations where an analytic likelihood approach is not possible.
<p>
Furthermore, as data and models grow in complexity, the two approaches can diverge greatly.
</section>
<section>
<h2 id="___sec19">Bayesian parameter estimation (multiple parameters, covariance) </h2>
<div class="alert alert-block alert-block alert-text-normal">
<b></b>
<ul>
<p><li> multidimensional posterior pdf:s</li>
<p><li> nuisance parameters (e.g. background subtraction?)</li>
<p><li> corner plots, covariance, correlations</li>
<p><li> best example?</li>
</ul>
</div>
</section>
<section>
<h2 id="___sec20">Bayesian model selection </h2>
<div class="alert alert-block alert-block alert-text-normal">
<b></b>
<ul>
<p><li> Bayesian evidence</li>
<p><li> Occam's razor</li>
<p><li> Best example? How many spectral lines are there?</li>
</ul>
</div>
</section>
</div> <!-- class="slides" -->
</div> <!-- class="reveal" -->
doc/pub/Bayesian/html/Bayesian-solarized.html
+492 -5
View File
@@ -68,7 +68,35 @@ div { text-align: justify; text-justify: inter-word; }
("Bayes' theorem", 2, None, '___sec5'),
("The friends of Bayes' theorem", 2, None, '___sec6'),
('Inference With Parametric Models', 2, None, '___sec7'),
('Illustrative examples with python code', 2, None, '___sec8')]}
('Illustrative examples with python code', 2, None, '___sec8'),
('Example: Is this a fair coin?', 2, None, '___sec9'),
('A few words on different priors', 2, None, '___sec10'),
('Bayesian parameter estimation (single parameter)',
2,
None,
'___sec11'),
('Example: Measured flux from a star', 3, None, '___sec12'),
('Simple Photon Counts: Frequentist Approach',
3,
None,
'___sec13'),
('Simple Photon Counts: Bayesian Approach', 3, None, '___sec14'),
('A note about priors', 3, None, '___sec15'),
('Simple Photon Counts: Bayesian approach in practice',
3,
None,
'___sec16'),
('Best estimates and confidence intervals', 3, None, '___sec17'),
('Simple Photon Counts: Best estimates and confidence intervals',
3,
None,
'___sec18'),
('Bayesian parameter estimation (multiple parameters, '
'covariance)',
2,
None,
'___sec19'),
('Bayesian model selection', 2, None, '___sec20')]}
end of tocinfo -->
<body>
@@ -115,7 +143,7 @@ MathJax.Hub.Config({
<center>[3] <b>Department of Physics and Astronomy and National Superconducting Cyclotron Laboratory, Michigan State University</b></center>
<br>
<p>
<center><h4>Mar 10, 2018</h4></center> <!-- date -->
<center><h4>Mar 13, 2018</h4></center> <!-- date -->
<br>
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
@@ -271,9 +299,10 @@ $$
The importance of this property to data analysis becomes apparent if we replace \( A \) and \( B \) by hypothesis(\( H \)) and data(\( D \)):
$$
$$
p(H|D,I) = \frac{p(D|H,I) p(H|I)}{p(D|I)}.
$$
\begin{align}
p(H|D,I) &= \frac{p(D|H,I) p(H|I)}{p(D|I)}.
\label{eq:bayes}
\end{align}
$$
The power of Bayes&#8217; theorem lies in the fact that it relates the quantity of interest, the probability that the hypothesis is true given the data, to the term we have a better chance of being able to assign, the probability that we would have observed the measured data if the hypothesis was true.
@@ -404,6 +433,464 @@ $$
</div>
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="___sec9">Example: Is this a fair coin? </h2>
Let us begin with the analysis of data from a simple coin-tossing experiment.
Given that we had observed 6 heads in 8 flips, would you think it was a fair coin? By fair, we mean that we would be prepared to lay an even 1 : 1 bet on the outcome of a flip being a head or a tail. If we decide that the coin was fair, the question which follows naturally is how sure are we that this was so; if it was not fair, how unfair do we think it was? Furthermore, if we were to continue collecting data for this particular coin, observing the outcomes of additional flips, how would we update our belief on the fairness of the coin?
<p>
A sensible way of formulating this problem is to consider a large number of hypotheses about the range in which the bias-weighting of the coin might lie. If we denote the bias-weighting by \( H \), then \( H = 0 \) and \( H = 1 \) can represent a coin which produces a tail or a head on every flip, respectively. There is a continuum of possibilities for the value of H between these limits, with \( H = 0.5 \) indicating a fair coin. Our state of knowledge about the fairness, or the degree of unfairness, of the coin is then completely summarized by specifying how much we believe these various propositions to be true.
<p>
Let us perform a computer simulation of a coin-tossing experiment. This provides the data that we will be analysing.
<p>
<!-- code=python (!bc pycod) typeset with pygments style "perldoc" -->
<div class="highlight" style="background: #eeeedd"><pre style="line-height: 125%"><span></span><span style="color: #8B008B; font-weight: bold">import</span> <span style="color: #008b45; text-decoration: underline">numpy</span> <span style="color: #8B008B; font-weight: bold">as</span> <span style="color: #008b45; text-decoration: underline">np</span>
<span style="color: #8B008B; font-weight: bold">import</span> <span style="color: #008b45; text-decoration: underline">matplotlib.pyplot</span> <span style="color: #8B008B; font-weight: bold">as</span> <span style="color: #008b45; text-decoration: underline">plt</span>
</pre></div>
<p>
<!-- code=python (!bc pycod) typeset with pygments style "perldoc" -->
<div class="highlight" style="background: #eeeedd"><pre style="line-height: 125%"><span></span>np.random.seed(<span style="color: #B452CD">999</span>) <span style="color: #228B22"># for reproducibility</span>
a=<span style="color: #B452CD">0.6</span> <span style="color: #228B22"># biased coin</span>
flips=np.random.rand(<span style="color: #B452CD">2</span>**<span style="color: #B452CD">12</span>) <span style="color: #228B22"># simulates 4096 coin flips</span>
heads=flips&lt;a <span style="color: #228B22"># boolean array, heads[i]=True if flip i is heads</span>
</pre></div>
<p>
In the light of this data, our inference about the fairness of this coin is summarized by the conditional pdf: \( p(H|D,I) \). This is, of course, shorthand for the limiting case of a continuum of propositions for the value of \( H \); that is to say, the probability that \( H \) lies in an infinitesimally narrow range is given by \( p(H|D,I) dH \).
<p>
To estimate this posterior pdf, we need to use Bayes&#8217; theorem \eqref{eq:bayes}. We will ignore the denominator \( p(D|I) \) as it does not involve bias-weighting explicitly, and it will therefore not affect the shape of the desired pdf. At the end we can evaluate the missing constant subsequently from the normalization condition
$$
\begin{equation}
\int_0^1 p(H|D,I) dH = 1.
\label{eq:coin_posterior_norm}
\end{equation}
$$
<p>
The prior pdf, \( p(H|I) \), represents what we know about the coin given only the information \( I \) that we are dealing with a &#8216;strange coin&#8217;. We could keep a very open mind about the nature of the coin; a simple probability assignment which reflects this is a uniform, or flat, prior
$$
\begin{equation}
p(H|I) = \left\{ \begin{array}{ll}
1 & 0 \le H \le 1, \\
0 & \mathrm{otherwise}.
\end{array} \right.
\label{eq:coin_prior_uniform}
\end{equation}
$$
We will get back later to the choice of prior and its effect on the analysis.
<p>
This prior state of knowledge, or ignorance, is modified by the data through the likelihood function \( p(D|H,I) \). It is a measure of the chance that we would have obtained the data that we actually observed, if the value of the bias-weighting was given (as known). If, in the conditioning information \( I \), we assume that the flips of the coin were independent events, so that the outcome of one did not influence that of another, then the probability of obtaining the data `R heads in N tosses' is given by the binomial distribution (we leave a formal definition of this to a statistics textbook)
$$
\begin{equation}
p(D|H,I) \propto H^R (1-H)^{N-R}.
\label{_auto1}
\end{equation}
$$
It seems reasonable because \( H \) is the chance of obtaining a head on any flip, and there were \( R \) of them, and \( 1-H \) is the corresponding probability for a tail, of which there were \( N-R \). We note that this binomial distribution also contains a normalization factor, but we will ignore it since it does not depend explicitly on \( H \), the quantity of interest. It will be absorbed by the normalization condition \eqref{eq:coin_posterior_norm}.
<p>
We perform the setup of this Bayesian framework on the computer.
<p>
<!-- code=python (!bc pycod) typeset with pygments style "perldoc" -->
<div class="highlight" style="background: #eeeedd"><pre style="line-height: 125%"><span></span><span style="color: #8B008B; font-weight: bold">def</span> <span style="color: #008b45">prior</span>(H):
p=np.zeros_like(H)
p[(<span style="color: #B452CD">0</span>&lt;=x)&amp;(x&lt;=<span style="color: #B452CD">1</span>)]=<span style="color: #B452CD">1</span> <span style="color: #228B22"># allowed range: 0&lt;=H&lt;=1</span>
<span style="color: #8B008B; font-weight: bold">return</span> p <span style="color: #228B22"># uniform prior</span>
<span style="color: #8B008B; font-weight: bold">def</span> <span style="color: #008b45">likelihood</span>(H,data):
N = <span style="color: #658b00">len</span>(data)
no_of_heads = <span style="color: #658b00">sum</span>(data)
no_of_tails = N - no_of_heads
<span style="color: #8B008B; font-weight: bold">return</span> H**no_of_heads * (<span style="color: #B452CD">1</span>-H)**no_of_tails
<span style="color: #8B008B; font-weight: bold">def</span> <span style="color: #008b45">posterior</span>(H,data):
p=prior(H)*likelihood(H,data)
norm=np.trapz(p,H)
<span style="color: #8B008B; font-weight: bold">return</span> p/norm
</pre></div>
<p>
The next step is to confront this setup with the simulated data. To get a feel for the result, it is instructive to see how the posterior pdf evolves as we obtain more and more data pertaining to the coin. The results of such an analyses is shown in Fig. <a href="#fig:coinflipping">1</a>.
<p>
<!-- code=python (!bc pycod) typeset with pygments style "perldoc" -->
<div class="highlight" style="background: #eeeedd"><pre style="line-height: 125%"><span></span>x=np.linspace(<span style="color: #B452CD">0</span>,<span style="color: #B452CD">1</span>,<span style="color: #B452CD">100</span>)
fig, axs = plt.subplots(nrows=<span style="color: #B452CD">4</span>,ncols=<span style="color: #B452CD">3</span>,sharex=<span style="color: #658b00">True</span>,sharey=<span style="color: #CD5555">&#39;row&#39;</span>)
axs_vec=np.reshape(axs,-<span style="color: #B452CD">1</span>)
axs_vec[<span style="color: #B452CD">0</span>].plot(x,prior(x))
<span style="color: #8B008B; font-weight: bold">for</span> ndouble <span style="color: #8B008B">in</span> <span style="color: #658b00">range</span>(<span style="color: #B452CD">11</span>):
ax=axs_vec[<span style="color: #B452CD">1</span>+ndouble]
ax.plot(x,posterior(x,heads[:<span style="color: #B452CD">2</span>**ndouble]))
ax.text(<span style="color: #B452CD">0.1</span>, <span style="color: #B452CD">0.8</span>, <span style="color: #CD5555">&#39;$N={0}$&#39;</span>.format(<span style="color: #B452CD">2</span>**ndouble), transform=ax.transAxes)
<span style="color: #8B008B; font-weight: bold">for</span> row <span style="color: #8B008B">in</span> <span style="color: #658b00">range</span>(<span style="color: #B452CD">4</span>): axs[row,<span style="color: #B452CD">0</span>].set_ylabel(<span style="color: #CD5555">&#39;$p(H|D_\mathrm{obs},I)$&#39;</span>)
<span style="color: #8B008B; font-weight: bold">for</span> col <span style="color: #8B008B">in</span> <span style="color: #658b00">range</span>(<span style="color: #B452CD">3</span>): axs[-<span style="color: #B452CD">1</span>,col].set_xlabel(<span style="color: #CD5555">&#39;$H$&#39;</span>)
</pre></div>
<p>
<center> <!-- figure label: --> <div id="fig:coinflipping"></div> <!-- FIGURE -->
<hr class="figure">
<center><p class="caption">Figure 1: The evolution of the posterior pdf for the bias-weighting of a coin, as the number of data available increases. The figure on the top left-hand corner of each panel shows the number of data included in the analysis. <!-- caption label: fig:coinflipping --> </p></center>
<p><img src="fig/coinflipping_fig_1.png" align="bottom" width=500></p>
</center>
<p>
The panel in the top left-hand corner shows the posterior pdf for \( H \) given no data, i.e., it is the same as the prior pdf of Eq. \eqref{eq:coin_prior_uniform}. It indicates that we have no more reason to believe that the coin is fair than we have to think that it is double-headed, double-tailed, or of any other intermediate bias-weighting.
<p>
The first flip is obviously tails. At this point we have no evidence that the coin has a side with heads, as indicated by the pdf going to zero as \( H \to 1 \). The second flip is obviously heads and we have now excluded both extreme options \( H=0 \) (double-tailed) and \( H=1 \) (double-headed). We can note that the posterior at this point has the simple form \( p(H|D,I) = H(1-H) \) for \( 0 \le H \le 1 \).
<p>
The remainder of Fig. <a href="#fig:coinflipping">1</a> shows how the posterior pdf evolves as the number of data analysed becomes larger and larger. We see that the position of the maximum moves around, but that the amount by which it does so decreases with the increasing number of observations. The width of the posterior pdf also becomes narrower with more data, indicating that we are becoming increasingly confident in our estimate of the bias-weighting. For the coin in this example, the best estimate of \( H \) eventually converges to 0.6, which, of course, was the value chosen to simulate the flips.
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="___sec10">A few words on different priors </h2>
<ul>
<li> uniform</li>
<li> Gaussian</li>
<li> Jeffrey's prior</li>
</ul>
Repeat the coin flipping experiment with other priors.
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="___sec11">Bayesian parameter estimation (single parameter) </h2>
<div class="alert alert-block alert-block alert-text-normal">
<b></b>
<p>
We will now consider the very important task of model parameter estimation using statistical inference.
<!-- begin inline comment -->
<font color="red">(<b>CF 1</b>: maybe stress that model parameters are not random variables, and the meaning of parameter estimation is therefore very different between frequentist and bayesian approaches.)</font>
<!-- end inline comment -->
<p>
Throughout this section we will consider a specific example that involves a model with a single parameter: &quot;Measured flux from a star&quot;.
</div>
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h3 id="___sec12">Example: Measured flux from a star </h3>
Adapted from the blog <a href="http://jakevdp.github.io" target="_blank">Pythonic Perambulations</a> by Jake VanderPlas.
<p>
Imagine that we point our telescope to the sky, and observe the light coming from a single star. For the time being, we'll assume that the star's true flux is constant with time, i.e. that is it has a fixed value \( F_\mathrm{true} \) (we'll also ignore effects like sky noise and other sources of systematic error). We'll assume that we perform a series of \( N \) measurements with our telescope, where the ith measurement reports the observed photon flux \( F_i \) and error \( e_i \) [<a id="link_footnote_2" href="#def_footnote_2">2</a>].
The question is, given this set of measurements \( D = \{F_i, e_i\} \), what is our best estimate of the true flux \( F_\mathrm{true} \)?
<p id="def_footnote_2"><a href="#link_footnote_2"><b>2:</b></a> We'll make the reasonable assumption that errors are Gaussian. In a Frequentist perspective, \( e_i \) is the standard deviation of the results of a single measurement event in the limit of repetitions of <em>that event</em>. In the Bayesian perspective, \( e_i \) is the standard deviation of the (Gaussian) probability distribution describing our knowledge of that particular measurement given its observed value.</p>
<p>
Because the measurements are number counts, a Poisson distribution is a good approximation to the measurement process:
<p>
<!-- code=python (!bc pycod) typeset with pygments style "perldoc" -->
<div class="highlight" style="background: #eeeedd"><pre style="line-height: 125%"><span></span>np.random.seed(<span style="color: #B452CD">1</span>) <span style="color: #228B22"># for repeatability</span>
F_true = <span style="color: #B452CD">1000</span> <span style="color: #228B22"># true flux, say number of photons measured in 1 second</span>
N = <span style="color: #B452CD">50</span> <span style="color: #228B22"># number of measurements</span>
F = stats.poisson(F_true).rvs(N)
<span style="color: #228B22"># N measurements of the flux</span>
e = np.sqrt(F) <span style="color: #228B22"># errors on Poisson counts estimated via square root</span>
</pre></div>
<p>
Now let's make a simple visualization of the &quot;observed&quot; data, see Fig. <a href="#fig:flux">2</a>.
<p>
<!-- code=python (!bc pycod) typeset with pygments style "perldoc" -->
<div class="highlight" style="background: #eeeedd"><pre style="line-height: 125%"><span></span>fig, ax = plt.subplots()
ax.errorbar(F, np.arange(N), xerr=e, fmt=<span style="color: #CD5555">&#39;ok&#39;</span>, ecolor=<span style="color: #CD5555">&#39;gray&#39;</span>, alpha=<span style="color: #B452CD">0.5</span>)
ax.vlines([F_true], <span style="color: #B452CD">0</span>, N, linewidth=<span style="color: #B452CD">5</span>, alpha=<span style="color: #B452CD">0.2</span>)
ax.set_xlabel(<span style="color: #CD5555">&quot;Flux&quot;</span>);ax.set_ylabel(<span style="color: #CD5555">&quot;measurement number&quot;</span>);
</pre></div>
<p>
<center> <!-- figure label: --> <div id="fig:flux"></div> <!-- FIGURE -->
<hr class="figure">
<center><p class="caption">Figure 2: Single photon counts (flux measurements). <!-- caption label: fig:flux --> </p></center>
<p><img src="fig/singlephotoncount_fig_1.png" align="bottom" width=400></p>
</center>
<p>
These measurements each have a different error \( e_i \) which is estimated from Poisson statistics using the standard square-root rule. In this toy example we already know the true flux \( F_\mathrm{true} \), but the question is this: given our measurements and errors, what is our best estimate of the true flux?
<p>
Let's take a look at the frequentist and Bayesian approaches to solving this.
<h3 id="___sec13">Simple Photon Counts: Frequentist Approach </h3>
We'll start with the classical frequentist maximum likelihood approach. Given a single observation \( D_i = (F_i, e_i) \), we can compute the probability distribution of the measurement given the true flux Ftrue given our assumption of Gaussian errors
$$
\begin{equation}
p(D_i | F_\mathrm{true}, I) = \frac{1}{\sqrt{2\pi e_i^2}} \exp \left( \frac{-(F_i-F_\mathrm{true})^2}{2e_i^2} \right).
\label{_auto2}
\end{equation}
$$
This should be read &quot;the probability of \( D_i \) given \( F_\mathrm{true} \)
equals ...&quot;. You should recognize this as a normal distribution with mean \( F_\mathrm{true} \) and standard deviation \( e_i \).
<p>
We construct the <em>likelihood function</em> by computing the product of the probabilities for each data point
$$
\begin{equation}
\mathcal{L}(D | F_\mathrm{true}, I) = \prod_{i=1}^N p(D_i | F_\mathrm{true}, I),
\label{_auto3}
\end{equation}
$$
here \( D = \{D_i\} \) represents the entire set of measurements. Because the value of the likelihood can become very small, it is often more convenient to instead compute the log-likelihood. Combining the previous two equations and computing the log, we have
$$
\begin{equation}
\log\mathcal{L} = -\frac{1}{2} \sum_{i=1}^N \left[ \log(2\pi e_i^2) + \frac{(F_i-F_\mathrm{true})^2}{e_i^2} \right].
\label{_auto4}
\end{equation}
$$
<p>
What we'd like to do is determine \( F_\mathrm{true} \) such that the likelihood is maximized. For this simple problem, the maximization can be computed analytically (i.e. by setting \( d\log\mathcal{L}/d F_\mathrm{true} = 0 \)). This results in the following observed estimate of \( F_\mathrm{true} \)
$$
\begin{equation}
F_\mathrm{est} = \sum_{i=1}^N w_i F_i; \quad w_i = 1/e_i^2.
\label{_auto5}
\end{equation}
$$
Notice that in the special case of all errors \( e_i \) being equal, this reduces to
$$
\begin{equation}
F_\mathrm{est} = \frac{1}{N} \sum_{i=1} F_i.
\label{_auto6}
\end{equation}
$$
That is, in agreement with intuition, \( F_\mathrm{est} \) is simply the mean of the observed data when errors are equal.
<p>
We can go further and ask what the error of our estimate is. In the frequentist approach, this can be accomplished by fitting a Gaussian approximation to the likelihood curve at maximum; in this simple case this can also be solved analytically (the sum of Gaussians is also a Gaussian). It can be shown that the standard deviation of this Gaussian approximation is
$$
\begin{equation}
\sigma_\mathrm{est} = \sum_{i=1}^N w_i.
\label{_auto7}
\end{equation}
$$
These results are fairly simple calculations; let's evaluate them for our toy dataset:
<p>
<!-- code=python (!bc pycod) typeset with pygments style "perldoc" -->
<div class="highlight" style="background: #eeeedd"><pre style="line-height: 125%"><span></span>w=<span style="color: #B452CD">1.</span>/e**<span style="color: #B452CD">2</span>
<span style="color: #8B008B; font-weight: bold">print</span>(<span style="color: #CD5555">&quot;&quot;&quot;</span>
<span style="color: #CD5555">F_true = {0}</span>
<span style="color: #CD5555">F_est = {1:.0f} +/- {2:.0f} (based on {3} measurements) &quot;&quot;&quot;</span>\
.format(F_true, (w * F).sum() / w.sum(), w.sum() ** -<span style="color: #B452CD">0.5</span>, N))
</pre></div>
<p>
<code>F_true = 1000</code> <br />
<code>F_est = 998 +/- 4 (based on 50 measurements)</code> <br />
<p>
We find that for 50 measurements of the flux, our estimate has an error of about 0.4% and is consistent with the input value.
<h3 id="___sec14">Simple Photon Counts: Bayesian Approach </h3>
The Bayesian approach, as you might expect, begins and ends with probabilities. Our hypothesis is that the star has a constant flux \( F_\mathrm{true} \). It recognizes that what we fundamentally want to compute is our knowledge of the parameters in question given the data and other information (such as our knowledge of uncertainties for the observed values), i.e. in this case, \( p(F_\mathrm{true} | D,I) \).
Note that this formulation of the problem is fundamentally contrary to the frequentist philosophy, which says that probabilities have no meaning for model parameters like \( F_\mathrm{true} \). Nevertheless, within the Bayesian philosophy this is perfectly acceptable.
<p>
To compute this result, Bayesians next apply Bayes' Theorem \eqref{eq:bayes}.
If we set the prior \( p(F_\mathrm{true}|I) \propto 1 \) (a flat prior), we find
\( p(F_\mathrm{true}|D,I) \propto p(D | F_\mathrm{true},I) \equiv \mathcal{L}(D | F_\mathrm{true},I) \)
and the Bayesian probability is maximized at precisely the same value as the frequentist result! So despite the philosophical differences, we see that (for this simple problem at least) the Bayesian and frequentist point estimates are equivalent.
<h3 id="___sec15">A note about priors </h3>
The prior allows inclusion of other information into the computation, which becomes very useful in cases where multiple measurement strategies are being combined to constrain a single model. The necessity to specify a prior, however, is one of the more controversial pieces of Bayesian analysis.
A frequentist will point out that the prior is problematic when no true prior information is available. Though it might seem straightforward to use a noninformative prior like the flat prior mentioned above, there are some <a href="http://normaldeviate.wordpress.com/2013/07/13/lost-causes-in-statistics-ii-noninformative- priors/comment-page-1/" target="_blank">surprisingly subtleties</a> involved. It turns out that in many situations, a truly noninformative prior does not exist! Frequentists point out that the subjective choice of a prior which necessarily biases your result has no place in statistical data analysis.
A Bayesian would counter that frequentism doesn't solve this problem, but simply skirts the question. Frequentism can often be viewed as simply a special case of the Bayesian approach for some (implicit) choice of the prior: a Bayesian would say that it's better to make this implicit choice explicit, even if the choice might include some subjectivity.
<h3 id="___sec16">Simple Photon Counts: Bayesian approach in practice </h3>
Leaving these philosophical debates aside for the time being, let's address how Bayesian results are generally computed in practice. For a one parameter problem like the one considered here, it's as simple as computing the posterior probability \( p(F_\mathrm{true} | D,I) \) as a function of \( F_\mathrm{true} \): this is the distribution reflecting our knowledge of the parameter \( F_\mathrm{true} \).
But as the dimension of the model grows, this direct approach becomes increasingly intractable. For this reason, Bayesian calculations often depend on sampling methods such as Markov Chain Monte Carlo (MCMC). For this practical example, let us apply an MCMC approach using Dan Foreman-Mackey's <a href="http://dan.iel.fm/emcee/current/" target="_blank">emcee</a> package. Keep in mind here that the goal is to generate a set of points drawn from the posterior probability distribution, and to use those points to determine the answer we seek.
To perform this MCMC, we start by defining Python functions for the prior \( p(F_\mathrm{true} | I) \), the likelihood \( p(D | F_\mathrm{true},I) \), and the posterior \( p(F_\mathrm{true} | D,I) \), noting that none of these need be properly normalized. Our model here is one-dimensional, but to handle multi-dimensional models we'll define the model in terms of an array of parameters \( \boldsymbol{\alpha} \), which in this case is \( \boldsymbol{\alpha} = [F_\mathrm{true}] \)
<p>
<!-- code=python (!bc pycod) typeset with pygments style "perldoc" -->
<div class="highlight" style="background: #eeeedd"><pre style="line-height: 125%"><span></span><span style="color: #8B008B; font-weight: bold">def</span> <span style="color: #008b45">log_prior</span>(alpha):
<span style="color: #8B008B; font-weight: bold">return</span> <span style="color: #B452CD">0</span> <span style="color: #228B22"># flat prior</span>
<span style="color: #8B008B; font-weight: bold">def</span> <span style="color: #008b45">log_likelihood</span>(alpha, F, e):
<span style="color: #8B008B; font-weight: bold">return</span> -<span style="color: #B452CD">0.5</span> * np.sum(np.log(<span style="color: #B452CD">2</span> * np.pi * e ** <span style="color: #B452CD">2</span>) \
+ (F - alpha[<span style="color: #B452CD">0</span>]) ** <span style="color: #B452CD">2</span> / e ** <span style="color: #B452CD">2</span>)
<span style="color: #8B008B; font-weight: bold">def</span> <span style="color: #008b45">log_posterior</span>(alpha, F, e):
<span style="color: #8B008B; font-weight: bold">return</span> log_prior(alpha) + log_likelihood(alpha, F, e)
</pre></div>
<p>
Now we set up the problem, including generating some random starting guesses for the multiple chains of points.
<p>
<!-- code=python (!bc pycod) typeset with pygments style "perldoc" -->
<div class="highlight" style="background: #eeeedd"><pre style="line-height: 125%"><span></span>ndim = <span style="color: #B452CD">1</span> <span style="color: #228B22"># number of parameters in the model</span>
nwalkers = <span style="color: #B452CD">50</span> <span style="color: #228B22"># number of MCMC walkers</span>
nburn = <span style="color: #B452CD">1000</span> <span style="color: #228B22"># &quot;burn-in&quot; period to let chains stabilize</span>
nsteps = <span style="color: #B452CD">2000</span> <span style="color: #228B22"># number of MCMC steps to take</span>
<span style="color: #228B22"># we&#39;ll start at random locations between 0 and 2000</span>
starting_guesses = <span style="color: #B452CD">2000</span> * np.random.rand(nwalkers, ndim)
sampler = emcee.EnsembleSampler(nwalkers, ndim, log_posterior, args=[F,e])
sampler.run_mcmc(starting_guesses, nsteps)
<span style="color: #228B22"># Shape of sampler.chain = (nwalkers, nsteps, ndim)</span>
<span style="color: #228B22"># Flatten the sampler chain and discard burn-in points:</span>
samples = sampler.chain[:, nburn:, :].reshape((-<span style="color: #B452CD">1</span>, ndim))
</pre></div>
<p>
If this all worked correctly, the array sample should contain a series of 50,000 points drawn from the posterior. Let's plot them and check. See results in Fig. <a href="#fig:flux-bayesian">3</a>.
<p>
<!-- code=python (!bc pycod) typeset with pygments style "perldoc" -->
<div class="highlight" style="background: #eeeedd"><pre style="line-height: 125%"><span></span>fig, ax = plt.subplots()
ax.hist(samples, bins=<span style="color: #B452CD">50</span>, histtype=<span style="color: #CD5555">&quot;stepfilled&quot;</span>, alpha=<span style="color: #B452CD">0.3</span>, normed=<span style="color: #658b00">True</span>)
ax.set_xlabel(<span style="color: #CD5555">r&#39;$F_\mathrm{est}$&#39;</span>)
ax.set_ylabel(<span style="color: #CD5555">r&#39;$p(F_\mathrm{est}|D,I)$&#39;</span>)
</pre></div>
<p>
<center> <!-- figure label: --> <div id="fig:flux-bayesian"></div> <!-- FIGURE -->
<hr class="figure">
<center><p class="caption">Figure 3: Bayesian posterior pdf (represented by a histogram of MCMC samples) from flux measurements. <!-- caption label: fig:flux-bayesian --> </p></center>
<p><img src="fig/singlephotoncount_fig_2.png" align="bottom" width=400></p>
</center>
<h3 id="___sec17">Best estimates and confidence intervals </h3>
The posterior distribution from our Bayesian data analysis is the key quantity that encodes our inference about the values of the model parameters, given the data and the relevant background information. Often, however, we wish to summarize this result with just a few numbers: the best estimate and a measure of its reliability.
<p>
There are a few different options for this. The choice of the most appropriate one depends mainly on the shape of the posterior distribution:
<p>
<em>Symmetric posterior pdfs</em>: Since the probability (density) associated with any particular value of the parameter is a measure of how much we believe that it lies in the neighbourhood of that point, our best estimate is given by the maximum of the posterior pdf. If we denote the quantity of interest by \( X \), with a posterior pdf \( P =p(X|D,I) \), then the best estimate of its value \( X_0 \) is given by the condition \( dP/dX|_{X=X_0}=0 \). Strictly speaking, we should also check the sign of the second derivative to ensure that \( X_0 \) represents a maximum.
<p>
To obtain a measure of the reliability of this best estimate, we need to look at the width or spread of the posterior pdf about \( X_0 \). When considering the behaviour of any function in the neighbourhood of a particular point, it is often helpful to carry out a Taylor series expansion; this is simply a standard tool for (locally) approximating a complicated function by a low-order polynomial. The linear term is zero at the maximum and the quadratic term is often the dominating one determining the width of the posterior pdf. Ignoring all the higher-order terms we arrive at the Gaussian approximation
$$
\begin{equation}
p(X|D,I) \approx \frac{1}{\sigma\sqrt{2\pi}} \exp \left[ -\frac{(x-\mu)^2}{2\sigma^2} \right],
\label{_auto8}
\end{equation}
$$
where the mean \( \mu = X_0 \) and the variance \( \sigma = \left( - \left. \frac{d^2L}{dX^2} \right|_{X_0} \right)^{-1/2} \), where \( L \) is the logarithm of the posterior \( P \). Our inference about the quantity of interest is conveyed very concisely, therefore, by the statement \( X = X_0 \pm \sigma \), and
$$
$$
p(X_0-\sigma < X < X_0+\sigma | D,I) = \int_{X_0-\sigma}^{X_0+\sigma} p(X|D,I) dX \approx 0.67.
$$
$$
<p>
<em>Asymmetric posterior pdfs</em>: While the maximum of the posterior (\( X_0 \)) can still be regarded as giving the best estimate, the true value is now more likely to be on one side of this rather than the other. Alternatively one can compute the mean value, \( \langle X \rangle = \int X p(X|D,I) dX \), although this tends to overemphasise very long tails. The best option is probably a compromise that can be employed when having access to a large sample from the posterior (as provided by an MCMC), namely to give the median of this ensamble.
<p>
Furthermore, the concept of an error-bar does not seem appropriate in this case, as it implicitly entails the idea of symmetry. A good way of expressing the reliability with which a parameter can be inferred, for an asymmetric posterior pdf, is rather through a <em>confidence interval</em>. Since the area under the posterior pdf between \( X_1 \) and \( X_2 \) is proportional to how much we believe that \( X \) lies in that range, the shortest interval that encloses 67% of the area represents a sensible measure of the uncertainty of the estimate. Obviously we can choose to provide some other degree-of-belief that we think is relevant for the case at hand. Assuming that the posterior pdf has been normalized, to have unit area, we need to find \( X_1 \) and \( X_2 \) such that:
$$
$$
p(X_1 < X < X_2 | D,I) = \int_{X_1}^{X_2} p(X|D,I) dX \approx 0.67,
$$
$$
where the difference \( X_2 - X_1 \) is as small as possible. The region \( X_1 < X < X_2 \) is then called the shortest 67% confidence interval.
<p>
<em>Multimodal posterior pdfs</em>: We can sometimes obtain posteriors which are multimodal; i.e. contains several disconnected regions with large probabilities. There is no difficulty when one of the maxima is very much larger than the others: we can simply ignore the subsidiary solutions, to a good approximation, and concentrate on the global maximum. The problem arises when there are several maxima of comparable magnitude. What do we now mean by a best estimate, and how should we quantify its reliability? The idea of a best estimate and an error-bar, or even a confidence interval, is merely an attempt to summarize the posterior with just two or three numbers; sometimes this just can&#8217;t be done, and so these concepts are not valid. For the bimodal case we might be able to characterize the posterior in terms of a few numbers: two best estimates and their associated error-bars, or disjoint confidence intervals. For a general multimodal pdf, the most honest thing we can do is just display the posterior itself.
<h3 id="___sec18">Simple Photon Counts: Best estimates and confidence intervals </h3>
To compute these numbers for our example, you would run:
<p>
<!-- code=python (!bc pycod) typeset with pygments style "perldoc" -->
<div class="highlight" style="background: #eeeedd"><pre style="line-height: 125%"><span></span>sampper=np.percentile(samples, [<span style="color: #B452CD">2.5</span>, <span style="color: #B452CD">16.5</span>, <span style="color: #B452CD">50</span>, <span style="color: #B452CD">83.5</span>, <span style="color: #B452CD">97.5</span>],axis=<span style="color: #B452CD">0</span>).flatten()
<span style="color: #8B008B; font-weight: bold">print</span>(<span style="color: #CD5555">&quot;&quot;&quot;</span>
<span style="color: #CD5555">F_true = {0}</span>
<span style="color: #CD5555">Based on {1} measurements the posterior point estimates are:</span>
<span style="color: #CD5555">...F_est = {2:.0f} +/- {3:.0f}</span>
<span style="color: #CD5555">or using credible intervals:</span>
<span style="color: #CD5555">...F_est = {4:.0f} (posterior median) </span>
<span style="color: #CD5555">...F_est in [{5:.0f}, {6:.0f}] (67% credible interval) </span>
<span style="color: #CD5555">...F_est in [{7:.0f}, {8:.0f}] (95% credible interval) &quot;&quot;&quot;</span>\
.format(F_true, N, np.mean(samples), np.std(samples), \
sampper[<span style="color: #B452CD">2</span>], sampper[<span style="color: #B452CD">1</span>], sampper[<span style="color: #B452CD">3</span>], sampper[<span style="color: #B452CD">0</span>], sampper[<span style="color: #B452CD">4</span>]))
</pre></div>
<p>
<code>F_true = 1000</code> <br />
<code>Based on 50 measurements the posterior point estimates are:</code> <br />
<code>...F_est = 998 +/- 4</code> <br />
<code>or using credible intervals:</code> <br />
<code>...F_est = 998 (posterior median)</code> <br />
<code>...F_est in [993, 1002] (67% credible interval)</code> <br />
<code>...F_est in [989, 1006] (95% credible interval)</code> <br />
<p>
In this particular example, the posterior pdf is actually a Gaussian (since it is constructed as a product of Gaussians), and the mean and variance from the quadratic approximation will agree exactly with the frequentist approach.
<p>
From this final result you might come away with the impression that the Bayesian method is unnecessarily complicated, and in this case it certainly is. Using an MCMC sampler to characterize a one-dimensional normal distribution is a bit like using the Death Star to destroy a beach ball, but we did this here because it demonstrates an approach that can scale to complicated posteriors in many, many dimensions, and can provide nice results in more complicated situations where an analytic likelihood approach is not possible.
<p>
Furthermore, as data and models grow in complexity, the two approaches can diverge greatly.
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="___sec19">Bayesian parameter estimation (multiple parameters, covariance) </h2>
<div class="alert alert-block alert-block alert-text-normal">
<b></b>
<p>
<ul>
<li> multidimensional posterior pdf:s</li>
<li> nuisance parameters (e.g. background subtraction?)</li>
<li> corner plots, covariance, correlations</li>
<li> best example?</li>
</ul>
</div>
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="___sec20">Bayesian model selection </h2>
<div class="alert alert-block alert-block alert-text-normal">
<b></b>
<p>
<ul>
<li> Bayesian evidence</li>
<li> Occam's razor</li>
<li> Best example? How many spectral lines are there?</li>
</ul>
</div>
<!-- ------------------- end of main content --------------- -->
doc/pub/Bayesian/html/Bayesian.html
+492 -5
View File
@@ -73,7 +73,35 @@ div { text-align: justify; text-justify: inter-word; }
("Bayes' theorem", 2, None, '___sec5'),
("The friends of Bayes' theorem", 2, None, '___sec6'),
('Inference With Parametric Models', 2, None, '___sec7'),
('Illustrative examples with python code', 2, None, '___sec8')]}
('Illustrative examples with python code', 2, None, '___sec8'),
('Example: Is this a fair coin?', 2, None, '___sec9'),
('A few words on different priors', 2, None, '___sec10'),
('Bayesian parameter estimation (single parameter)',
2,
None,
'___sec11'),
('Example: Measured flux from a star', 3, None, '___sec12'),
('Simple Photon Counts: Frequentist Approach',
3,
None,
'___sec13'),
('Simple Photon Counts: Bayesian Approach', 3, None, '___sec14'),
('A note about priors', 3, None, '___sec15'),
('Simple Photon Counts: Bayesian approach in practice',
3,
None,
'___sec16'),
('Best estimates and confidence intervals', 3, None, '___sec17'),
('Simple Photon Counts: Best estimates and confidence intervals',
3,
None,
'___sec18'),
('Bayesian parameter estimation (multiple parameters, '
'covariance)',
2,
None,
'___sec19'),
('Bayesian model selection', 2, None, '___sec20')]}
end of tocinfo -->
<body>
@@ -120,7 +148,7 @@ MathJax.Hub.Config({
<center>[3] <b>Department of Physics and Astronomy and National Superconducting Cyclotron Laboratory, Michigan State University</b></center>
<br>
<p>
<center><h4>Mar 10, 2018</h4></center> <!-- date -->
<center><h4>Mar 13, 2018</h4></center> <!-- date -->
<br>
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
@@ -276,9 +304,10 @@ $$
The importance of this property to data analysis becomes apparent if we replace \( A \) and \( B \) by hypothesis(\( H \)) and data(\( D \)):
$$
$$
p(H|D,I) = \frac{p(D|H,I) p(H|I)}{p(D|I)}.
$$
\begin{align}
p(H|D,I) &= \frac{p(D|H,I) p(H|I)}{p(D|I)}.
\label{eq:bayes}
\end{align}
$$
The power of Bayes&#8217; theorem lies in the fact that it relates the quantity of interest, the probability that the hypothesis is true given the data, to the term we have a better chance of being able to assign, the probability that we would have observed the measured data if the hypothesis was true.
@@ -409,6 +438,464 @@ $$
</div>
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="___sec9">Example: Is this a fair coin? </h2>
Let us begin with the analysis of data from a simple coin-tossing experiment.
Given that we had observed 6 heads in 8 flips, would you think it was a fair coin? By fair, we mean that we would be prepared to lay an even 1 : 1 bet on the outcome of a flip being a head or a tail. If we decide that the coin was fair, the question which follows naturally is how sure are we that this was so; if it was not fair, how unfair do we think it was? Furthermore, if we were to continue collecting data for this particular coin, observing the outcomes of additional flips, how would we update our belief on the fairness of the coin?
<p>
A sensible way of formulating this problem is to consider a large number of hypotheses about the range in which the bias-weighting of the coin might lie. If we denote the bias-weighting by \( H \), then \( H = 0 \) and \( H = 1 \) can represent a coin which produces a tail or a head on every flip, respectively. There is a continuum of possibilities for the value of H between these limits, with \( H = 0.5 \) indicating a fair coin. Our state of knowledge about the fairness, or the degree of unfairness, of the coin is then completely summarized by specifying how much we believe these various propositions to be true.
<p>
Let us perform a computer simulation of a coin-tossing experiment. This provides the data that we will be analysing.
<p>
<!-- code=python (!bc pycod) typeset with pygments style "default" -->
<div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%"><span></span><span style="color: #008000; font-weight: bold">import</span> <span style="color: #0000FF; font-weight: bold">numpy</span> <span style="color: #008000; font-weight: bold">as</span> <span style="color: #0000FF; font-weight: bold">np</span>
<span style="color: #008000; font-weight: bold">import</span> <span style="color: #0000FF; font-weight: bold">matplotlib.pyplot</span> <span style="color: #008000; font-weight: bold">as</span> <span style="color: #0000FF; font-weight: bold">plt</span>
</pre></div>
<p>
<!-- code=python (!bc pycod) typeset with pygments style "default" -->
<div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%"><span></span>np<span style="color: #666666">.</span>random<span style="color: #666666">.</span>seed(<span style="color: #666666">999</span>) <span style="color: #408080; font-style: italic"># for reproducibility</span>
a<span style="color: #666666">=0.6</span> <span style="color: #408080; font-style: italic"># biased coin</span>
flips<span style="color: #666666">=</span>np<span style="color: #666666">.</span>random<span style="color: #666666">.</span>rand(<span style="color: #666666">2**12</span>) <span style="color: #408080; font-style: italic"># simulates 4096 coin flips</span>
heads<span style="color: #666666">=</span>flips<span style="color: #666666">&lt;</span>a <span style="color: #408080; font-style: italic"># boolean array, heads[i]=True if flip i is heads</span>
</pre></div>
<p>
In the light of this data, our inference about the fairness of this coin is summarized by the conditional pdf: \( p(H|D,I) \). This is, of course, shorthand for the limiting case of a continuum of propositions for the value of \( H \); that is to say, the probability that \( H \) lies in an infinitesimally narrow range is given by \( p(H|D,I) dH \).
<p>
To estimate this posterior pdf, we need to use Bayes&#8217; theorem \eqref{eq:bayes}. We will ignore the denominator \( p(D|I) \) as it does not involve bias-weighting explicitly, and it will therefore not affect the shape of the desired pdf. At the end we can evaluate the missing constant subsequently from the normalization condition
$$
\begin{equation}
\int_0^1 p(H|D,I) dH = 1.
\label{eq:coin_posterior_norm}
\end{equation}
$$
<p>
The prior pdf, \( p(H|I) \), represents what we know about the coin given only the information \( I \) that we are dealing with a &#8216;strange coin&#8217;. We could keep a very open mind about the nature of the coin; a simple probability assignment which reflects this is a uniform, or flat, prior
$$
\begin{equation}
p(H|I) = \left\{ \begin{array}{ll}
1 & 0 \le H \le 1, \\
0 & \mathrm{otherwise}.
\end{array} \right.
\label{eq:coin_prior_uniform}
\end{equation}
$$
We will get back later to the choice of prior and its effect on the analysis.
<p>
This prior state of knowledge, or ignorance, is modified by the data through the likelihood function \( p(D|H,I) \). It is a measure of the chance that we would have obtained the data that we actually observed, if the value of the bias-weighting was given (as known). If, in the conditioning information \( I \), we assume that the flips of the coin were independent events, so that the outcome of one did not influence that of another, then the probability of obtaining the data `R heads in N tosses' is given by the binomial distribution (we leave a formal definition of this to a statistics textbook)
$$
\begin{equation}
p(D|H,I) \propto H^R (1-H)^{N-R}.
\label{_auto1}
\end{equation}
$$
It seems reasonable because \( H \) is the chance of obtaining a head on any flip, and there were \( R \) of them, and \( 1-H \) is the corresponding probability for a tail, of which there were \( N-R \). We note that this binomial distribution also contains a normalization factor, but we will ignore it since it does not depend explicitly on \( H \), the quantity of interest. It will be absorbed by the normalization condition \eqref{eq:coin_posterior_norm}.
<p>
We perform the setup of this Bayesian framework on the computer.
<p>
<!-- code=python (!bc pycod) typeset with pygments style "default" -->
<div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%"><span></span><span style="color: #008000; font-weight: bold">def</span> <span style="color: #0000FF">prior</span>(H):
p<span style="color: #666666">=</span>np<span style="color: #666666">.</span>zeros_like(H)
p[(<span style="color: #666666">0&lt;=</span>x)<span style="color: #666666">&amp;</span>(x<span style="color: #666666">&lt;=1</span>)]<span style="color: #666666">=1</span> <span style="color: #408080; font-style: italic"># allowed range: 0&lt;=H&lt;=1</span>
<span style="color: #008000; font-weight: bold">return</span> p <span style="color: #408080; font-style: italic"># uniform prior</span>
<span style="color: #008000; font-weight: bold">def</span> <span style="color: #0000FF">likelihood</span>(H,data):
N <span style="color: #666666">=</span> <span style="color: #008000">len</span>(data)
no_of_heads <span style="color: #666666">=</span> <span style="color: #008000">sum</span>(data)
no_of_tails <span style="color: #666666">=</span> N <span style="color: #666666">-</span> no_of_heads
<span style="color: #008000; font-weight: bold">return</span> H<span style="color: #666666">**</span>no_of_heads <span style="color: #666666">*</span> (<span style="color: #666666">1-</span>H)<span style="color: #666666">**</span>no_of_tails
<span style="color: #008000; font-weight: bold">def</span> <span style="color: #0000FF">posterior</span>(H,data):
p<span style="color: #666666">=</span>prior(H)<span style="color: #666666">*</span>likelihood(H,data)
norm<span style="color: #666666">=</span>np<span style="color: #666666">.</span>trapz(p,H)
<span style="color: #008000; font-weight: bold">return</span> p<span style="color: #666666">/</span>norm
</pre></div>
<p>
The next step is to confront this setup with the simulated data. To get a feel for the result, it is instructive to see how the posterior pdf evolves as we obtain more and more data pertaining to the coin. The results of such an analyses is shown in Fig. <a href="#fig:coinflipping">1</a>.
<p>
<!-- code=python (!bc pycod) typeset with pygments style "default" -->
<div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%"><span></span>x<span style="color: #666666">=</span>np<span style="color: #666666">.</span>linspace(<span style="color: #666666">0</span>,<span style="color: #666666">1</span>,<span style="color: #666666">100</span>)
fig, axs <span style="color: #666666">=</span> plt<span style="color: #666666">.</span>subplots(nrows<span style="color: #666666">=4</span>,ncols<span style="color: #666666">=3</span>,sharex<span style="color: #666666">=</span><span style="color: #008000">True</span>,sharey<span style="color: #666666">=</span><span style="color: #BA2121">&#39;row&#39;</span>)
axs_vec<span style="color: #666666">=</span>np<span style="color: #666666">.</span>reshape(axs,<span style="color: #666666">-1</span>)
axs_vec[<span style="color: #666666">0</span>]<span style="color: #666666">.</span>plot(x,prior(x))
<span style="color: #008000; font-weight: bold">for</span> ndouble <span style="color: #AA22FF; font-weight: bold">in</span> <span style="color: #008000">range</span>(<span style="color: #666666">11</span>):
ax<span style="color: #666666">=</span>axs_vec[<span style="color: #666666">1+</span>ndouble]
ax<span style="color: #666666">.</span>plot(x,posterior(x,heads[:<span style="color: #666666">2**</span>ndouble]))
ax<span style="color: #666666">.</span>text(<span style="color: #666666">0.1</span>, <span style="color: #666666">0.8</span>, <span style="color: #BA2121">&#39;$N={0}$&#39;</span><span style="color: #666666">.</span>format(<span style="color: #666666">2**</span>ndouble), transform<span style="color: #666666">=</span>ax<span style="color: #666666">.</span>transAxes)
<span style="color: #008000; font-weight: bold">for</span> row <span style="color: #AA22FF; font-weight: bold">in</span> <span style="color: #008000">range</span>(<span style="color: #666666">4</span>): axs[row,<span style="color: #666666">0</span>]<span style="color: #666666">.</span>set_ylabel(<span style="color: #BA2121">&#39;$p(H|D_\mathrm{obs},I)$&#39;</span>)
<span style="color: #008000; font-weight: bold">for</span> col <span style="color: #AA22FF; font-weight: bold">in</span> <span style="color: #008000">range</span>(<span style="color: #666666">3</span>): axs[<span style="color: #666666">-1</span>,col]<span style="color: #666666">.</span>set_xlabel(<span style="color: #BA2121">&#39;$H$&#39;</span>)
</pre></div>
<p>
<center> <!-- figure label: --> <div id="fig:coinflipping"></div> <!-- FIGURE -->
<hr class="figure">
<center><p class="caption">Figure 1: The evolution of the posterior pdf for the bias-weighting of a coin, as the number of data available increases. The figure on the top left-hand corner of each panel shows the number of data included in the analysis. <!-- caption label: fig:coinflipping --> </p></center>
<p><img src="fig/coinflipping_fig_1.png" align="bottom" width=500></p>
</center>
<p>
The panel in the top left-hand corner shows the posterior pdf for \( H \) given no data, i.e., it is the same as the prior pdf of Eq. \eqref{eq:coin_prior_uniform}. It indicates that we have no more reason to believe that the coin is fair than we have to think that it is double-headed, double-tailed, or of any other intermediate bias-weighting.
<p>
The first flip is obviously tails. At this point we have no evidence that the coin has a side with heads, as indicated by the pdf going to zero as \( H \to 1 \). The second flip is obviously heads and we have now excluded both extreme options \( H=0 \) (double-tailed) and \( H=1 \) (double-headed). We can note that the posterior at this point has the simple form \( p(H|D,I) = H(1-H) \) for \( 0 \le H \le 1 \).
<p>
The remainder of Fig. <a href="#fig:coinflipping">1</a> shows how the posterior pdf evolves as the number of data analysed becomes larger and larger. We see that the position of the maximum moves around, but that the amount by which it does so decreases with the increasing number of observations. The width of the posterior pdf also becomes narrower with more data, indicating that we are becoming increasingly confident in our estimate of the bias-weighting. For the coin in this example, the best estimate of \( H \) eventually converges to 0.6, which, of course, was the value chosen to simulate the flips.
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="___sec10">A few words on different priors </h2>
<ul>
<li> uniform</li>
<li> Gaussian</li>
<li> Jeffrey's prior</li>
</ul>
Repeat the coin flipping experiment with other priors.
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="___sec11">Bayesian parameter estimation (single parameter) </h2>
<div class="alert alert-block alert-block alert-text-normal">
<b></b>
<p>
We will now consider the very important task of model parameter estimation using statistical inference.
<!-- begin inline comment -->
<font color="red">(<b>CF 1</b>: maybe stress that model parameters are not random variables, and the meaning of parameter estimation is therefore very different between frequentist and bayesian approaches.)</font>
<!-- end inline comment -->
<p>
Throughout this section we will consider a specific example that involves a model with a single parameter: &quot;Measured flux from a star&quot;.
</div>
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h3 id="___sec12">Example: Measured flux from a star </h3>
Adapted from the blog <a href="http://jakevdp.github.io" target="_blank">Pythonic Perambulations</a> by Jake VanderPlas.
<p>
Imagine that we point our telescope to the sky, and observe the light coming from a single star. For the time being, we'll assume that the star's true flux is constant with time, i.e. that is it has a fixed value \( F_\mathrm{true} \) (we'll also ignore effects like sky noise and other sources of systematic error). We'll assume that we perform a series of \( N \) measurements with our telescope, where the ith measurement reports the observed photon flux \( F_i \) and error \( e_i \) [<a id="link_footnote_2" href="#def_footnote_2">2</a>].
The question is, given this set of measurements \( D = \{F_i, e_i\} \), what is our best estimate of the true flux \( F_\mathrm{true} \)?
<p id="def_footnote_2"><a href="#link_footnote_2"><b>2:</b></a> We'll make the reasonable assumption that errors are Gaussian. In a Frequentist perspective, \( e_i \) is the standard deviation of the results of a single measurement event in the limit of repetitions of <em>that event</em>. In the Bayesian perspective, \( e_i \) is the standard deviation of the (Gaussian) probability distribution describing our knowledge of that particular measurement given its observed value.</p>
<p>
Because the measurements are number counts, a Poisson distribution is a good approximation to the measurement process:
<p>
<!-- code=python (!bc pycod) typeset with pygments style "default" -->
<div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%"><span></span>np<span style="color: #666666">.</span>random<span style="color: #666666">.</span>seed(<span style="color: #666666">1</span>) <span style="color: #408080; font-style: italic"># for repeatability</span>
F_true <span style="color: #666666">=</span> <span style="color: #666666">1000</span> <span style="color: #408080; font-style: italic"># true flux, say number of photons measured in 1 second</span>
N <span style="color: #666666">=</span> <span style="color: #666666">50</span> <span style="color: #408080; font-style: italic"># number of measurements</span>
F <span style="color: #666666">=</span> stats<span style="color: #666666">.</span>poisson(F_true)<span style="color: #666666">.</span>rvs(N)
<span style="color: #408080; font-style: italic"># N measurements of the flux</span>
e <span style="color: #666666">=</span> np<span style="color: #666666">.</span>sqrt(F) <span style="color: #408080; font-style: italic"># errors on Poisson counts estimated via square root</span>
</pre></div>
<p>
Now let's make a simple visualization of the &quot;observed&quot; data, see Fig. <a href="#fig:flux">2</a>.
<p>
<!-- code=python (!bc pycod) typeset with pygments style "default" -->
<div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%"><span></span>fig, ax <span style="color: #666666">=</span> plt<span style="color: #666666">.</span>subplots()
ax<span style="color: #666666">.</span>errorbar(F, np<span style="color: #666666">.</span>arange(N), xerr<span style="color: #666666">=</span>e, fmt<span style="color: #666666">=</span><span style="color: #BA2121">&#39;ok&#39;</span>, ecolor<span style="color: #666666">=</span><span style="color: #BA2121">&#39;gray&#39;</span>, alpha<span style="color: #666666">=0.5</span>)
ax<span style="color: #666666">.</span>vlines([F_true], <span style="color: #666666">0</span>, N, linewidth<span style="color: #666666">=5</span>, alpha<span style="color: #666666">=0.2</span>)
ax<span style="color: #666666">.</span>set_xlabel(<span style="color: #BA2121">&quot;Flux&quot;</span>);ax<span style="color: #666666">.</span>set_ylabel(<span style="color: #BA2121">&quot;measurement number&quot;</span>);
</pre></div>
<p>
<center> <!-- figure label: --> <div id="fig:flux"></div> <!-- FIGURE -->
<hr class="figure">
<center><p class="caption">Figure 2: Single photon counts (flux measurements). <!-- caption label: fig:flux --> </p></center>
<p><img src="fig/singlephotoncount_fig_1.png" align="bottom" width=400></p>
</center>
<p>
These measurements each have a different error \( e_i \) which is estimated from Poisson statistics using the standard square-root rule. In this toy example we already know the true flux \( F_\mathrm{true} \), but the question is this: given our measurements and errors, what is our best estimate of the true flux?
<p>
Let's take a look at the frequentist and Bayesian approaches to solving this.
<h3 id="___sec13">Simple Photon Counts: Frequentist Approach </h3>
We'll start with the classical frequentist maximum likelihood approach. Given a single observation \( D_i = (F_i, e_i) \), we can compute the probability distribution of the measurement given the true flux Ftrue given our assumption of Gaussian errors
$$
\begin{equation}
p(D_i | F_\mathrm{true}, I) = \frac{1}{\sqrt{2\pi e_i^2}} \exp \left( \frac{-(F_i-F_\mathrm{true})^2}{2e_i^2} \right).
\label{_auto2}
\end{equation}
$$
This should be read &quot;the probability of \( D_i \) given \( F_\mathrm{true} \)
equals ...&quot;. You should recognize this as a normal distribution with mean \( F_\mathrm{true} \) and standard deviation \( e_i \).
<p>
We construct the <em>likelihood function</em> by computing the product of the probabilities for each data point
$$
\begin{equation}
\mathcal{L}(D | F_\mathrm{true}, I) = \prod_{i=1}^N p(D_i | F_\mathrm{true}, I),
\label{_auto3}
\end{equation}
$$
here \( D = \{D_i\} \) represents the entire set of measurements. Because the value of the likelihood can become very small, it is often more convenient to instead compute the log-likelihood. Combining the previous two equations and computing the log, we have
$$
\begin{equation}
\log\mathcal{L} = -\frac{1}{2} \sum_{i=1}^N \left[ \log(2\pi e_i^2) + \frac{(F_i-F_\mathrm{true})^2}{e_i^2} \right].
\label{_auto4}
\end{equation}
$$
<p>
What we'd like to do is determine \( F_\mathrm{true} \) such that the likelihood is maximized. For this simple problem, the maximization can be computed analytically (i.e. by setting \( d\log\mathcal{L}/d F_\mathrm{true} = 0 \)). This results in the following observed estimate of \( F_\mathrm{true} \)
$$
\begin{equation}
F_\mathrm{est} = \sum_{i=1}^N w_i F_i; \quad w_i = 1/e_i^2.
\label{_auto5}
\end{equation}
$$
Notice that in the special case of all errors \( e_i \) being equal, this reduces to
$$
\begin{equation}
F_\mathrm{est} = \frac{1}{N} \sum_{i=1} F_i.
\label{_auto6}
\end{equation}
$$
That is, in agreement with intuition, \( F_\mathrm{est} \) is simply the mean of the observed data when errors are equal.
<p>
We can go further and ask what the error of our estimate is. In the frequentist approach, this can be accomplished by fitting a Gaussian approximation to the likelihood curve at maximum; in this simple case this can also be solved analytically (the sum of Gaussians is also a Gaussian). It can be shown that the standard deviation of this Gaussian approximation is
$$
\begin{equation}
\sigma_\mathrm{est} = \sum_{i=1}^N w_i.
\label{_auto7}
\end{equation}
$$
These results are fairly simple calculations; let's evaluate them for our toy dataset:
<p>
<!-- code=python (!bc pycod) typeset with pygments style "default" -->
<div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%"><span></span>w<span style="color: #666666">=1./</span>e<span style="color: #666666">**2</span>
<span style="color: #008000; font-weight: bold">print</span>(<span style="color: #BA2121">&quot;&quot;&quot;</span>
<span style="color: #BA2121">F_true = {0}</span>
<span style="color: #BA2121">F_est = {1:.0f} +/- {2:.0f} (based on {3} measurements) &quot;&quot;&quot;</span>\
<span style="color: #666666">.</span>format(F_true, (w <span style="color: #666666">*</span> F)<span style="color: #666666">.</span>sum() <span style="color: #666666">/</span> w<span style="color: #666666">.</span>sum(), w<span style="color: #666666">.</span>sum() <span style="color: #666666">**</span> <span style="color: #666666">-0.5</span>, N))
</pre></div>
<p>
<code>F_true = 1000</code> <br />
<code>F_est = 998 +/- 4 (based on 50 measurements)</code> <br />
<p>
We find that for 50 measurements of the flux, our estimate has an error of about 0.4% and is consistent with the input value.
<h3 id="___sec14">Simple Photon Counts: Bayesian Approach </h3>
The Bayesian approach, as you might expect, begins and ends with probabilities. Our hypothesis is that the star has a constant flux \( F_\mathrm{true} \). It recognizes that what we fundamentally want to compute is our knowledge of the parameters in question given the data and other information (such as our knowledge of uncertainties for the observed values), i.e. in this case, \( p(F_\mathrm{true} | D,I) \).
Note that this formulation of the problem is fundamentally contrary to the frequentist philosophy, which says that probabilities have no meaning for model parameters like \( F_\mathrm{true} \). Nevertheless, within the Bayesian philosophy this is perfectly acceptable.
<p>
To compute this result, Bayesians next apply Bayes' Theorem \eqref{eq:bayes}.
If we set the prior \( p(F_\mathrm{true}|I) \propto 1 \) (a flat prior), we find
\( p(F_\mathrm{true}|D,I) \propto p(D | F_\mathrm{true},I) \equiv \mathcal{L}(D | F_\mathrm{true},I) \)
and the Bayesian probability is maximized at precisely the same value as the frequentist result! So despite the philosophical differences, we see that (for this simple problem at least) the Bayesian and frequentist point estimates are equivalent.
<h3 id="___sec15">A note about priors </h3>
The prior allows inclusion of other information into the computation, which becomes very useful in cases where multiple measurement strategies are being combined to constrain a single model. The necessity to specify a prior, however, is one of the more controversial pieces of Bayesian analysis.
A frequentist will point out that the prior is problematic when no true prior information is available. Though it might seem straightforward to use a noninformative prior like the flat prior mentioned above, there are some <a href="http://normaldeviate.wordpress.com/2013/07/13/lost-causes-in-statistics-ii-noninformative- priors/comment-page-1/" target="_blank">surprisingly subtleties</a> involved. It turns out that in many situations, a truly noninformative prior does not exist! Frequentists point out that the subjective choice of a prior which necessarily biases your result has no place in statistical data analysis.
A Bayesian would counter that frequentism doesn't solve this problem, but simply skirts the question. Frequentism can often be viewed as simply a special case of the Bayesian approach for some (implicit) choice of the prior: a Bayesian would say that it's better to make this implicit choice explicit, even if the choice might include some subjectivity.
<h3 id="___sec16">Simple Photon Counts: Bayesian approach in practice </h3>
Leaving these philosophical debates aside for the time being, let's address how Bayesian results are generally computed in practice. For a one parameter problem like the one considered here, it's as simple as computing the posterior probability \( p(F_\mathrm{true} | D,I) \) as a function of \( F_\mathrm{true} \): this is the distribution reflecting our knowledge of the parameter \( F_\mathrm{true} \).
But as the dimension of the model grows, this direct approach becomes increasingly intractable. For this reason, Bayesian calculations often depend on sampling methods such as Markov Chain Monte Carlo (MCMC). For this practical example, let us apply an MCMC approach using Dan Foreman-Mackey's <a href="http://dan.iel.fm/emcee/current/" target="_blank">emcee</a> package. Keep in mind here that the goal is to generate a set of points drawn from the posterior probability distribution, and to use those points to determine the answer we seek.
To perform this MCMC, we start by defining Python functions for the prior \( p(F_\mathrm{true} | I) \), the likelihood \( p(D | F_\mathrm{true},I) \), and the posterior \( p(F_\mathrm{true} | D,I) \), noting that none of these need be properly normalized. Our model here is one-dimensional, but to handle multi-dimensional models we'll define the model in terms of an array of parameters \( \boldsymbol{\alpha} \), which in this case is \( \boldsymbol{\alpha} = [F_\mathrm{true}] \)
<p>
<!-- code=python (!bc pycod) typeset with pygments style "default" -->
<div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%"><span></span><span style="color: #008000; font-weight: bold">def</span> <span style="color: #0000FF">log_prior</span>(alpha):
<span style="color: #008000; font-weight: bold">return</span> <span style="color: #666666">0</span> <span style="color: #408080; font-style: italic"># flat prior</span>
<span style="color: #008000; font-weight: bold">def</span> <span style="color: #0000FF">log_likelihood</span>(alpha, F, e):
<span style="color: #008000; font-weight: bold">return</span> <span style="color: #666666">-0.5</span> <span style="color: #666666">*</span> np<span style="color: #666666">.</span>sum(np<span style="color: #666666">.</span>log(<span style="color: #666666">2</span> <span style="color: #666666">*</span> np<span style="color: #666666">.</span>pi <span style="color: #666666">*</span> e <span style="color: #666666">**</span> <span style="color: #666666">2</span>) \
<span style="color: #666666">+</span> (F <span style="color: #666666">-</span> alpha[<span style="color: #666666">0</span>]) <span style="color: #666666">**</span> <span style="color: #666666">2</span> <span style="color: #666666">/</span> e <span style="color: #666666">**</span> <span style="color: #666666">2</span>)
<span style="color: #008000; font-weight: bold">def</span> <span style="color: #0000FF">log_posterior</span>(alpha, F, e):
<span style="color: #008000; font-weight: bold">return</span> log_prior(alpha) <span style="color: #666666">+</span> log_likelihood(alpha, F, e)
</pre></div>
<p>
Now we set up the problem, including generating some random starting guesses for the multiple chains of points.
<p>
<!-- code=python (!bc pycod) typeset with pygments style "default" -->
<div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%"><span></span>ndim <span style="color: #666666">=</span> <span style="color: #666666">1</span> <span style="color: #408080; font-style: italic"># number of parameters in the model</span>
nwalkers <span style="color: #666666">=</span> <span style="color: #666666">50</span> <span style="color: #408080; font-style: italic"># number of MCMC walkers</span>
nburn <span style="color: #666666">=</span> <span style="color: #666666">1000</span> <span style="color: #408080; font-style: italic"># &quot;burn-in&quot; period to let chains stabilize</span>
nsteps <span style="color: #666666">=</span> <span style="color: #666666">2000</span> <span style="color: #408080; font-style: italic"># number of MCMC steps to take</span>
<span style="color: #408080; font-style: italic"># we&#39;ll start at random locations between 0 and 2000</span>
starting_guesses <span style="color: #666666">=</span> <span style="color: #666666">2000</span> <span style="color: #666666">*</span> np<span style="color: #666666">.</span>random<span style="color: #666666">.</span>rand(nwalkers, ndim)
sampler <span style="color: #666666">=</span> emcee<span style="color: #666666">.</span>EnsembleSampler(nwalkers, ndim, log_posterior, args<span style="color: #666666">=</span>[F,e])
sampler<span style="color: #666666">.</span>run_mcmc(starting_guesses, nsteps)
<span style="color: #408080; font-style: italic"># Shape of sampler.chain = (nwalkers, nsteps, ndim)</span>
<span style="color: #408080; font-style: italic"># Flatten the sampler chain and discard burn-in points:</span>
samples <span style="color: #666666">=</span> sampler<span style="color: #666666">.</span>chain[:, nburn:, :]<span style="color: #666666">.</span>reshape((<span style="color: #666666">-1</span>, ndim))
</pre></div>
<p>
If this all worked correctly, the array sample should contain a series of 50,000 points drawn from the posterior. Let's plot them and check. See results in Fig. <a href="#fig:flux-bayesian">3</a>.
<p>
<!-- code=python (!bc pycod) typeset with pygments style "default" -->
<div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%"><span></span>fig, ax <span style="color: #666666">=</span> plt<span style="color: #666666">.</span>subplots()
ax<span style="color: #666666">.</span>hist(samples, bins<span style="color: #666666">=50</span>, histtype<span style="color: #666666">=</span><span style="color: #BA2121">&quot;stepfilled&quot;</span>, alpha<span style="color: #666666">=0.3</span>, normed<span style="color: #666666">=</span><span style="color: #008000">True</span>)
ax<span style="color: #666666">.</span>set_xlabel(<span style="color: #BA2121">r&#39;$F_\mathrm{est}$&#39;</span>)
ax<span style="color: #666666">.</span>set_ylabel(<span style="color: #BA2121">r&#39;$p(F_\mathrm{est}|D,I)$&#39;</span>)
</pre></div>
<p>
<center> <!-- figure label: --> <div id="fig:flux-bayesian"></div> <!-- FIGURE -->
<hr class="figure">
<center><p class="caption">Figure 3: Bayesian posterior pdf (represented by a histogram of MCMC samples) from flux measurements. <!-- caption label: fig:flux-bayesian --> </p></center>
<p><img src="fig/singlephotoncount_fig_2.png" align="bottom" width=400></p>
</center>
<h3 id="___sec17">Best estimates and confidence intervals </h3>
The posterior distribution from our Bayesian data analysis is the key quantity that encodes our inference about the values of the model parameters, given the data and the relevant background information. Often, however, we wish to summarize this result with just a few numbers: the best estimate and a measure of its reliability.
<p>
There are a few different options for this. The choice of the most appropriate one depends mainly on the shape of the posterior distribution:
<p>
<em>Symmetric posterior pdfs</em>: Since the probability (density) associated with any particular value of the parameter is a measure of how much we believe that it lies in the neighbourhood of that point, our best estimate is given by the maximum of the posterior pdf. If we denote the quantity of interest by \( X \), with a posterior pdf \( P =p(X|D,I) \), then the best estimate of its value \( X_0 \) is given by the condition \( dP/dX|_{X=X_0}=0 \). Strictly speaking, we should also check the sign of the second derivative to ensure that \( X_0 \) represents a maximum.
<p>
To obtain a measure of the reliability of this best estimate, we need to look at the width or spread of the posterior pdf about \( X_0 \). When considering the behaviour of any function in the neighbourhood of a particular point, it is often helpful to carry out a Taylor series expansion; this is simply a standard tool for (locally) approximating a complicated function by a low-order polynomial. The linear term is zero at the maximum and the quadratic term is often the dominating one determining the width of the posterior pdf. Ignoring all the higher-order terms we arrive at the Gaussian approximation
$$
\begin{equation}
p(X|D,I) \approx \frac{1}{\sigma\sqrt{2\pi}} \exp \left[ -\frac{(x-\mu)^2}{2\sigma^2} \right],
\label{_auto8}
\end{equation}
$$
where the mean \( \mu = X_0 \) and the variance \( \sigma = \left( - \left. \frac{d^2L}{dX^2} \right|_{X_0} \right)^{-1/2} \), where \( L \) is the logarithm of the posterior \( P \). Our inference about the quantity of interest is conveyed very concisely, therefore, by the statement \( X = X_0 \pm \sigma \), and
$$
$$
p(X_0-\sigma < X < X_0+\sigma | D,I) = \int_{X_0-\sigma}^{X_0+\sigma} p(X|D,I) dX \approx 0.67.
$$
$$
<p>
<em>Asymmetric posterior pdfs</em>: While the maximum of the posterior (\( X_0 \)) can still be regarded as giving the best estimate, the true value is now more likely to be on one side of this rather than the other. Alternatively one can compute the mean value, \( \langle X \rangle = \int X p(X|D,I) dX \), although this tends to overemphasise very long tails. The best option is probably a compromise that can be employed when having access to a large sample from the posterior (as provided by an MCMC), namely to give the median of this ensamble.
<p>
Furthermore, the concept of an error-bar does not seem appropriate in this case, as it implicitly entails the idea of symmetry. A good way of expressing the reliability with which a parameter can be inferred, for an asymmetric posterior pdf, is rather through a <em>confidence interval</em>. Since the area under the posterior pdf between \( X_1 \) and \( X_2 \) is proportional to how much we believe that \( X \) lies in that range, the shortest interval that encloses 67% of the area represents a sensible measure of the uncertainty of the estimate. Obviously we can choose to provide some other degree-of-belief that we think is relevant for the case at hand. Assuming that the posterior pdf has been normalized, to have unit area, we need to find \( X_1 \) and \( X_2 \) such that:
$$
$$
p(X_1 < X < X_2 | D,I) = \int_{X_1}^{X_2} p(X|D,I) dX \approx 0.67,
$$
$$
where the difference \( X_2 - X_1 \) is as small as possible. The region \( X_1 < X < X_2 \) is then called the shortest 67% confidence interval.
<p>
<em>Multimodal posterior pdfs</em>: We can sometimes obtain posteriors which are multimodal; i.e. contains several disconnected regions with large probabilities. There is no difficulty when one of the maxima is very much larger than the others: we can simply ignore the subsidiary solutions, to a good approximation, and concentrate on the global maximum. The problem arises when there are several maxima of comparable magnitude. What do we now mean by a best estimate, and how should we quantify its reliability? The idea of a best estimate and an error-bar, or even a confidence interval, is merely an attempt to summarize the posterior with just two or three numbers; sometimes this just can&#8217;t be done, and so these concepts are not valid. For the bimodal case we might be able to characterize the posterior in terms of a few numbers: two best estimates and their associated error-bars, or disjoint confidence intervals. For a general multimodal pdf, the most honest thing we can do is just display the posterior itself.
<h3 id="___sec18">Simple Photon Counts: Best estimates and confidence intervals </h3>
To compute these numbers for our example, you would run:
<p>
<!-- code=python (!bc pycod) typeset with pygments style "default" -->
<div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%"><span></span>sampper<span style="color: #666666">=</span>np<span style="color: #666666">.</span>percentile(samples, [<span style="color: #666666">2.5</span>, <span style="color: #666666">16.5</span>, <span style="color: #666666">50</span>, <span style="color: #666666">83.5</span>, <span style="color: #666666">97.5</span>],axis<span style="color: #666666">=0</span>)<span style="color: #666666">.</span>flatten()
<span style="color: #008000; font-weight: bold">print</span>(<span style="color: #BA2121">&quot;&quot;&quot;</span>
<span style="color: #BA2121">F_true = {0}</span>
<span style="color: #BA2121">Based on {1} measurements the posterior point estimates are:</span>
<span style="color: #BA2121">...F_est = {2:.0f} +/- {3:.0f}</span>
<span style="color: #BA2121">or using credible intervals:</span>
<span style="color: #BA2121">...F_est = {4:.0f} (posterior median) </span>
<span style="color: #BA2121">...F_est in [{5:.0f}, {6:.0f}] (67</span><span style="color: #BB6688; font-weight: bold">% c</span><span style="color: #BA2121">redible interval) </span>
<span style="color: #BA2121">...F_est in [{7:.0f}, {8:.0f}] (95</span><span style="color: #BB6688; font-weight: bold">% c</span><span style="color: #BA2121">redible interval) &quot;&quot;&quot;</span>\
<span style="color: #666666">.</span>format(F_true, N, np<span style="color: #666666">.</span>mean(samples), np<span style="color: #666666">.</span>std(samples), \
sampper[<span style="color: #666666">2</span>], sampper[<span style="color: #666666">1</span>], sampper[<span style="color: #666666">3</span>], sampper[<span style="color: #666666">0</span>], sampper[<span style="color: #666666">4</span>]))
</pre></div>
<p>
<code>F_true = 1000</code> <br />
<code>Based on 50 measurements the posterior point estimates are:</code> <br />
<code>...F_est = 998 +/- 4</code> <br />
<code>or using credible intervals:</code> <br />
<code>...F_est = 998 (posterior median)</code> <br />
<code>...F_est in [993, 1002] (67% credible interval)</code> <br />
<code>...F_est in [989, 1006] (95% credible interval)</code> <br />
<p>
In this particular example, the posterior pdf is actually a Gaussian (since it is constructed as a product of Gaussians), and the mean and variance from the quadratic approximation will agree exactly with the frequentist approach.
<p>
From this final result you might come away with the impression that the Bayesian method is unnecessarily complicated, and in this case it certainly is. Using an MCMC sampler to characterize a one-dimensional normal distribution is a bit like using the Death Star to destroy a beach ball, but we did this here because it demonstrates an approach that can scale to complicated posteriors in many, many dimensions, and can provide nice results in more complicated situations where an analytic likelihood approach is not possible.
<p>
Furthermore, as data and models grow in complexity, the two approaches can diverge greatly.
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="___sec19">Bayesian parameter estimation (multiple parameters, covariance) </h2>
<div class="alert alert-block alert-block alert-text-normal">
<b></b>
<p>
<ul>
<li> multidimensional posterior pdf:s</li>
<li> nuisance parameters (e.g. background subtraction?)</li>
<li> corner plots, covariance, correlations</li>
<li> best example?</li>
</ul>
</div>
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="___sec20">Bayesian model selection </h2>
<div class="alert alert-block alert-block alert-text-normal">
<b></b>
<p>
<ul>
<li> Bayesian evidence</li>
<li> Occam's razor</li>
<li> Best example? How many spectral lines are there?</li>
</ul>
</div>
<!-- ------------------- end of main content --------------- -->
doc/pub/Bayesian/ipynb/Bayesian.ipynb
+706 -5
View File
@@ -12,7 +12,7 @@
"<!-- dom:AUTHOR: Morten Hjorth-Jensen at Department of Physics, University of Oslo & Department of Physics and Astronomy and National Superconducting Cyclotron Laboratory, Michigan State University -->\n",
"<!-- Author: --> **Morten Hjorth-Jensen**, Department of Physics, University of Oslo and Department of Physics and Astronomy and National Superconducting Cyclotron Laboratory, Michigan State University\n",
"\n",
"Date: **Mar 10, 2018**\n",
"Date: **Mar 13, 2018**\n",
"\n",
"Copyright 1999-2018, Morten Hjorth-Jensen. Released under CC Attribution-NonCommercial 4.0 license\n",
"\n",
@@ -158,10 +158,14 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"eq:bayes\"></div>\n",
"\n",
"$$\n",
"$$\n",
"\\begin{equation}\n",
"p(H|D,I) = \\frac{p(D|H,I) p(H|I)}{p(D|I)}.\n",
"$$\n",
"\\label{eq:bayes} \\tag{1}\n",
"\\end{equation}\n",
"$$"
]
},
@@ -282,11 +286,708 @@
"\n",
"* Linear fit with outliers (nuisance parameters)\n",
"\n",
"* ..."
"* ...\n",
"\n",
"\n",
"\n",
"## Example: Is this a fair coin?\n",
"Let us begin with the analysis of data from a simple coin-tossing experiment. \n",
"Given that we had observed 6 heads in 8 flips, would you think it was a fair coin? By fair, we mean that we would be prepared to lay an even 1 : 1 bet on the outcome of a flip being a head or a tail. If we decide that the coin was fair, the question which follows naturally is how sure are we that this was so; if it was not fair, how unfair do we think it was? Furthermore, if we were to continue collecting data for this particular coin, observing the outcomes of additional flips, how would we update our belief on the fairness of the coin?\n",
"\n",
"A sensible way of formulating this problem is to consider a large number of hypotheses about the range in which the bias-weighting of the coin might lie. If we denote the bias-weighting by $H$, then $H = 0$ and $H = 1$ can represent a coin which produces a tail or a head on every flip, respectively. There is a continuum of possibilities for the value of H between these limits, with $H = 0.5$ indicating a fair coin. Our state of knowledge about the fairness, or the degree of unfairness, of the coin is then completely summarized by specifying how much we believe these various propositions to be true. \n",
"\n",
"Let us perform a computer simulation of a coin-tossing experiment. This provides the data that we will be analysing."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"0\n",
" \n",
"<\n",
"<\n",
"<\n",
"!\n",
"!\n",
"C\n",
"O\n",
"D\n",
"E\n",
"_\n",
"B\n",
"L\n",
"O\n",
"C\n",
"K\n",
" \n",
" \n",
"p\n",
"y\n",
"c\n",
"o\n",
"d"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"np.random.seed(999) # for reproducibility\n",
"a=0.6 # biased coin\n",
"flips=np.random.rand(2**12) # simulates 4096 coin flips\n",
"heads=flips<a # boolean array, heads[i]=True if flip i is heads"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In the light of this data, our inference about the fairness of this coin is summarized by the conditional pdf: $p(H|D,I)$. This is, of course, shorthand for the limiting case of a continuum of propositions for the value of $H$; that is to say, the probability that $H$ lies in an infinitesimally narrow range is given by $p(H|D,I) dH$. \n",
"\n",
"To estimate this posterior pdf, we need to use Bayes’ theorem ([eq:bayes](#eq:bayes)). We will ignore the denominator $p(D|I)$ as it does not involve bias-weighting explicitly, and it will therefore not affect the shape of the desired pdf. At the end we can evaluate the missing constant subsequently from the normalization condition"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"eq:coin_posterior_norm\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"\\int_0^1 p(H|D,I) dH = 1.\n",
"\\label{eq:coin_posterior_norm} \\tag{2}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The prior pdf, $p(H|I)$, represents what we know about the coin given only the information $I$ that we are dealing with a ‘strange coin’. We could keep a very open mind about the nature of the coin; a simple probability assignment which reflects this is a uniform, or flat, prior"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"eq:coin_prior_uniform\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"p(H|I) = \\left\\{ \\begin{array}{ll}\n",
"1 & 0 \\le H \\le 1, \\\\\n",
"0 & \\mathrm{otherwise}.\n",
"\\end{array} \\right.\n",
"\\label{eq:coin_prior_uniform} \\tag{3}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We will get back later to the choice of prior and its effect on the analysis.\n",
"\n",
"This prior state of knowledge, or ignorance, is modified by the data through the likelihood function $p(D|H,I)$. It is a measure of the chance that we would have obtained the data that we actually observed, if the value of the bias-weighting was given (as known). If, in the conditioning information $I$, we assume that the flips of the coin were independent events, so that the outcome of one did not influence that of another, then the probability of obtaining the data `R heads in N tosses' is given by the binomial distribution (we leave a formal definition of this to a statistics textbook)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"_auto1\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"p(D|H,I) \\propto H^R (1-H)^{N-R}.\n",
"\\label{_auto1} \\tag{4}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"It seems reasonable because $H$ is the chance of obtaining a head on any flip, and there were $R$ of them, and $1-H$ is the corresponding probability for a tail, of which there were $N-R$. We note that this binomial distribution also contains a normalization factor, but we will ignore it since it does not depend explicitly on $H$, the quantity of interest. It will be absorbed by the normalization condition ([eq:coin_posterior_norm](#eq:coin_posterior_norm)).\n",
"\n",
"We perform the setup of this Bayesian framework on the computer."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"def prior(H):\n",
" p=np.zeros_like(H)\n",
" p[(0<=x)&(x<=1)]=1 # allowed range: 0<=H<=1\n",
" return p # uniform prior\n",
"def likelihood(H,data):\n",
" N = len(data)\n",
" no_of_heads = sum(data)\n",
" no_of_tails = N - no_of_heads\n",
" return H**no_of_heads * (1-H)**no_of_tails\n",
"def posterior(H,data):\n",
" p=prior(H)*likelihood(H,data)\n",
" norm=np.trapz(p,H)\n",
" return p/norm"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The next step is to confront this setup with the simulated data. To get a feel for the result, it is instructive to see how the posterior pdf evolves as we obtain more and more data pertaining to the coin. The results of such an analyses is shown in Fig. [fig:coinflipping](#fig:coinflipping)."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"x=np.linspace(0,1,100)\n",
"fig, axs = plt.subplots(nrows=4,ncols=3,sharex=True,sharey='row')\n",
"axs_vec=np.reshape(axs,-1)\n",
"axs_vec[0].plot(x,prior(x))\n",
"for ndouble in range(11):\n",
" ax=axs_vec[1+ndouble]\n",
" ax.plot(x,posterior(x,heads[:2**ndouble]))\n",
" ax.text(0.1, 0.8, '$N={0}$'.format(2**ndouble), transform=ax.transAxes)\n",
"for row in range(4): axs[row,0].set_ylabel('$p(H|D_\\mathrm{obs},I)$')\n",
"for col in range(3): axs[-1,col].set_xlabel('$H$')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- dom:FIGURE:[fig/coinflipping_fig_1.png, width=500 frac=0.95] The evolution of the posterior pdf for the bias-weighting of a coin, as the number of data available increases. The figure on the top left-hand corner of each panel shows the number of data included in the analysis. <div id=\"fig:coinflipping\"></div> -->\n",
"<!-- begin figure -->\n",
"<div id=\"fig:coinflipping\"></div>\n",
"\n",
"<p>The evolution of the posterior pdf for the bias-weighting of a coin, as the number of data available increases. The figure on the top left-hand corner of each panel shows the number of data included in the analysis.</p>\n",
"<img src=\"fig/coinflipping_fig_1.png\" width=500>\n",
"\n",
"<!-- end figure -->\n",
"\n",
"\n",
"The panel in the top left-hand corner shows the posterior pdf for $H$ given no data, i.e., it is the same as the prior pdf of Eq. ([eq:coin_prior_uniform](#eq:coin_prior_uniform)). It indicates that we have no more reason to believe that the coin is fair than we have to think that it is double-headed, double-tailed, or of any other intermediate bias-weighting.\n",
"\n",
"The first flip is obviously tails. At this point we have no evidence that the coin has a side with heads, as indicated by the pdf going to zero as $H \\to 1$. The second flip is obviously heads and we have now excluded both extreme options $H=0$ (double-tailed) and $H=1$ (double-headed). We can note that the posterior at this point has the simple form $p(H|D,I) = H(1-H)$ for $0 \\le H \\le 1$.\n",
"\n",
"The remainder of Fig. [fig:coinflipping](#fig:coinflipping) shows how the posterior pdf evolves as the number of data analysed becomes larger and larger. We see that the position of the maximum moves around, but that the amount by which it does so decreases with the increasing number of observations. The width of the posterior pdf also becomes narrower with more data, indicating that we are becoming increasingly confident in our estimate of the bias-weighting. For the coin in this example, the best estimate of $H$ eventually converges to 0.6, which, of course, was the value chosen to simulate the flips.\n",
"\n",
"## A few words on different priors\n",
"* uniform\n",
"\n",
"* Gaussian\n",
"\n",
"* Jeffrey's prior\n",
"\n",
"Repeat the coin flipping experiment with other priors.\n",
"\n",
"## Bayesian parameter estimation (single parameter)\n",
"We will now consider the very important task of model parameter estimation using statistical inference. \n",
"[CF 1: maybe stress that model parameters are not random variables, and the meaning of parameter estimation is therefore very different between frequentist and bayesian approaches.]\n",
"\n",
"Throughout this section we will consider a specific example that involves a model with a single parameter: \"Measured flux from a star\".\n",
"\n",
"\n",
"\n",
"### Example: Measured flux from a star\n",
"\n",
"Adapted from the blog [Pythonic Perambulations](http://jakevdp.github.io) by Jake VanderPlas.\n",
"\n",
"Imagine that we point our telescope to the sky, and observe the light coming from a single star. For the time being, we'll assume that the star's true flux is constant with time, i.e. that is it has a fixed value $F_\\mathrm{true}$ (we'll also ignore effects like sky noise and other sources of systematic error). We'll assume that we perform a series of $N$ measurements with our telescope, where the ith measurement reports the observed photon flux $F_i$ and error $e_i$[^errors].\n",
"The question is, given this set of measurements $D = \\{F_i, e_i\\}$, what is our best estimate of the true flux $F_\\mathrm{true}$?\n",
"\n",
"[^errors]: We'll make the reasonable assumption that errors are Gaussian. In a Frequentist perspective, $e_i$ is the standard deviation of the results of a single measurement event in the limit of repetitions of *that event*. In the Bayesian perspective, $e_i$ is the standard deviation of the (Gaussian) probability distribution describing our knowledge of that particular measurement given its observed value.\n",
"\n",
"Because the measurements are number counts, a Poisson distribution is a good approximation to the measurement process:"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"np.random.seed(1) # for repeatability\n",
"F_true = 1000 # true flux, say number of photons measured in 1 second\n",
"N = 50 # number of measurements\n",
"F = stats.poisson(F_true).rvs(N)\n",
" # N measurements of the flux\n",
"e = np.sqrt(F) # errors on Poisson counts estimated via square root"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now let's make a simple visualization of the \"observed\" data, see Fig. [fig:flux](#fig:flux)."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"fig, ax = plt.subplots()\n",
"ax.errorbar(F, np.arange(N), xerr=e, fmt='ok', ecolor='gray', alpha=0.5)\n",
"ax.vlines([F_true], 0, N, linewidth=5, alpha=0.2)\n",
"ax.set_xlabel(\"Flux\");ax.set_ylabel(\"measurement number\");"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- dom:FIGURE:[fig/singlephotoncount_fig_1.png, width=400 frac=0.8] Single photon counts (flux measurements). <div id=\"fig:flux\"></div> -->\n",
"<!-- begin figure -->\n",
"<div id=\"fig:flux\"></div>\n",
"\n",
"<p>Single photon counts (flux measurements).</p>\n",
"<img src=\"fig/singlephotoncount_fig_1.png\" width=400>\n",
"\n",
"<!-- end figure -->\n",
"\n",
"\n",
"These measurements each have a different error $e_i$ which is estimated from Poisson statistics using the standard square-root rule. In this toy example we already know the true flux $F_\\mathrm{true}$, but the question is this: given our measurements and errors, what is our best estimate of the true flux?\n",
"\n",
"Let's take a look at the frequentist and Bayesian approaches to solving this.\n",
"\n",
"### Simple Photon Counts: Frequentist Approach\n",
"\n",
"We'll start with the classical frequentist maximum likelihood approach. Given a single observation $D_i = (F_i, e_i)$, we can compute the probability distribution of the measurement given the true flux Ftrue given our assumption of Gaussian errors"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"_auto2\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"p(D_i | F_\\mathrm{true}, I) = \\frac{1}{\\sqrt{2\\pi e_i^2}} \\exp \\left( \\frac{-(F_i-F_\\mathrm{true})^2}{2e_i^2} \\right).\n",
"\\label{_auto2} \\tag{5}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This should be read \"the probability of $D_i$ given $F_\\mathrm{true}$\n",
"equals ...\". You should recognize this as a normal distribution with mean $F_\\mathrm{true}$ and standard deviation $e_i$.\n",
"\n",
"We construct the *likelihood function* by computing the product of the probabilities for each data point"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"_auto3\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"\\mathcal{L}(D | F_\\mathrm{true}, I) = \\prod_{i=1}^N p(D_i | F_\\mathrm{true}, I),\n",
"\\label{_auto3} \\tag{6}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"here $D = \\{D_i\\}$ represents the entire set of measurements. Because the value of the likelihood can become very small, it is often more convenient to instead compute the log-likelihood. Combining the previous two equations and computing the log, we have"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"_auto4\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"\\log\\mathcal{L} = -\\frac{1}{2} \\sum_{i=1}^N \\left[ \\log(2\\pi e_i^2) + \\frac{(F_i-F_\\mathrm{true})^2}{e_i^2} \\right].\n",
"\\label{_auto4} \\tag{7}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"What we'd like to do is determine $F_\\mathrm{true}$ such that the likelihood is maximized. For this simple problem, the maximization can be computed analytically (i.e. by setting $d\\log\\mathcal{L}/d F_\\mathrm{true} = 0$). This results in the following observed estimate of $F_\\mathrm{true}$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"_auto5\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"F_\\mathrm{est} = \\sum_{i=1}^N w_i F_i; \\quad w_i = 1/e_i^2.\n",
"\\label{_auto5} \\tag{8}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Notice that in the special case of all errors $e_i$ being equal, this reduces to"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"_auto6\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"F_\\mathrm{est} = \\frac{1}{N} \\sum_{i=1} F_i.\n",
"\\label{_auto6} \\tag{9}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"That is, in agreement with intuition, $F_\\mathrm{est}$ is simply the mean of the observed data when errors are equal.\n",
"\n",
"We can go further and ask what the error of our estimate is. In the frequentist approach, this can be accomplished by fitting a Gaussian approximation to the likelihood curve at maximum; in this simple case this can also be solved analytically (the sum of Gaussians is also a Gaussian). It can be shown that the standard deviation of this Gaussian approximation is"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"_auto7\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"\\sigma_\\mathrm{est} = \\sum_{i=1}^N w_i.\n",
"\\label{_auto7} \\tag{10}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"These results are fairly simple calculations; let's evaluate them for our toy dataset:"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"w=1./e**2\n",
"print(\"\"\"\n",
"F_true = {0}\n",
"F_est = {1:.0f} +/- {2:.0f} (based on {3} measurements) \"\"\"\\\n",
" .format(F_true, (w * F).sum() / w.sum(), w.sum() ** -0.5, N))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"`F_true = 1000` \n",
"`F_est = 998 +/- 4 (based on 50 measurements)` \n",
"\n",
"We find that for 50 measurements of the flux, our estimate has an error of about 0.4% and is consistent with the input value.\n",
"\n",
"\n",
"### Simple Photon Counts: Bayesian Approach\n",
"\n",
"The Bayesian approach, as you might expect, begins and ends with probabilities. Our hypothesis is that the star has a constant flux $F_\\mathrm{true}$. It recognizes that what we fundamentally want to compute is our knowledge of the parameters in question given the data and other information (such as our knowledge of uncertainties for the observed values), i.e. in this case, $p(F_\\mathrm{true} | D,I)$.\n",
"Note that this formulation of the problem is fundamentally contrary to the frequentist philosophy, which says that probabilities have no meaning for model parameters like $F_\\mathrm{true}$. Nevertheless, within the Bayesian philosophy this is perfectly acceptable.\n",
"\n",
"To compute this result, Bayesians next apply Bayes' Theorem ([eq:bayes](#eq:bayes)).\n",
"If we set the prior $p(F_\\mathrm{true}|I) \\propto 1$ (a flat prior), we find\n",
"$p(F_\\mathrm{true}|D,I) \\propto p(D | F_\\mathrm{true},I) \\equiv \\mathcal{L}(D | F_\\mathrm{true},I)$\n",
"and the Bayesian probability is maximized at precisely the same value as the frequentist result! So despite the philosophical differences, we see that (for this simple problem at least) the Bayesian and frequentist point estimates are equivalent.\n",
"\n",
"### A note about priors\n",
"\n",
"The prior allows inclusion of other information into the computation, which becomes very useful in cases where multiple measurement strategies are being combined to constrain a single model. The necessity to specify a prior, however, is one of the more controversial pieces of Bayesian analysis.\n",
"A frequentist will point out that the prior is problematic when no true prior information is available. Though it might seem straightforward to use a noninformative prior like the flat prior mentioned above, there are some [surprisingly subtleties](http://normaldeviate.wordpress.com/2013/07/13/lost-causes-in-statistics-ii-noninformative- priors/comment-page-1/) involved. It turns out that in many situations, a truly noninformative prior does not exist! Frequentists point out that the subjective choice of a prior which necessarily biases your result has no place in statistical data analysis.\n",
"A Bayesian would counter that frequentism doesn't solve this problem, but simply skirts the question. Frequentism can often be viewed as simply a special case of the Bayesian approach for some (implicit) choice of the prior: a Bayesian would say that it's better to make this implicit choice explicit, even if the choice might include some subjectivity.\n",
"\n",
"### Simple Photon Counts: Bayesian approach in practice\n",
"\n",
"Leaving these philosophical debates aside for the time being, let's address how Bayesian results are generally computed in practice. For a one parameter problem like the one considered here, it's as simple as computing the posterior probability $p(F_\\mathrm{true} | D,I)$ as a function of $F_\\mathrm{true}$: this is the distribution reflecting our knowledge of the parameter $F_\\mathrm{true}$.\n",
"But as the dimension of the model grows, this direct approach becomes increasingly intractable. For this reason, Bayesian calculations often depend on sampling methods such as Markov Chain Monte Carlo (MCMC). For this practical example, let us apply an MCMC approach using Dan Foreman-Mackey's [emcee](http://dan.iel.fm/emcee/current/) package. Keep in mind here that the goal is to generate a set of points drawn from the posterior probability distribution, and to use those points to determine the answer we seek.\n",
"To perform this MCMC, we start by defining Python functions for the prior $p(F_\\mathrm{true} | I)$, the likelihood $p(D | F_\\mathrm{true},I)$, and the posterior $p(F_\\mathrm{true} | D,I)$, noting that none of these need be properly normalized. Our model here is one-dimensional, but to handle multi-dimensional models we'll define the model in terms of an array of parameters $\\boldsymbol{\\alpha}$, which in this case is $\\boldsymbol{\\alpha} = [F_\\mathrm{true}]$"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"def log_prior(alpha):\n",
" return 0 # flat prior\n",
"\n",
"def log_likelihood(alpha, F, e):\n",
" return -0.5 * np.sum(np.log(2 * np.pi * e ** 2) \\\n",
" + (F - alpha[0]) ** 2 / e ** 2)\n",
" \n",
"def log_posterior(alpha, F, e):\n",
" return log_prior(alpha) + log_likelihood(alpha, F, e)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we set up the problem, including generating some random starting guesses for the multiple chains of points."
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"ndim = 1 # number of parameters in the model\n",
"nwalkers = 50 # number of MCMC walkers\n",
"nburn = 1000 # \"burn-in\" period to let chains stabilize\n",
"nsteps = 2000 # number of MCMC steps to take\n",
"# we'll start at random locations between 0 and 2000\n",
"starting_guesses = 2000 * np.random.rand(nwalkers, ndim)\n",
"sampler = emcee.EnsembleSampler(nwalkers, ndim, log_posterior, args=[F,e])\n",
"sampler.run_mcmc(starting_guesses, nsteps)\n",
"# Shape of sampler.chain = (nwalkers, nsteps, ndim)\n",
"# Flatten the sampler chain and discard burn-in points:\n",
"samples = sampler.chain[:, nburn:, :].reshape((-1, ndim))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"If this all worked correctly, the array sample should contain a series of 50,000 points drawn from the posterior. Let's plot them and check. See results in Fig. [fig:flux-bayesian](#fig:flux-bayesian)."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"fig, ax = plt.subplots()\n",
"ax.hist(samples, bins=50, histtype=\"stepfilled\", alpha=0.3, normed=True)\n",
"ax.set_xlabel(r'$F_\\mathrm{est}$')\n",
"ax.set_ylabel(r'$p(F_\\mathrm{est}|D,I)$')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- dom:FIGURE:[fig/singlephotoncount_fig_2.png, width=400 frac=0.8] Bayesian posterior pdf (represented by a histogram of MCMC samples) from flux measurements. <div id=\"fig:flux-bayesian\"></div> -->\n",
"<!-- begin figure -->\n",
"<div id=\"fig:flux-bayesian\"></div>\n",
"\n",
"<p>Bayesian posterior pdf (represented by a histogram of MCMC samples) from flux measurements.</p>\n",
"<img src=\"fig/singlephotoncount_fig_2.png\" width=400>\n",
"\n",
"<!-- end figure -->\n",
"\n",
"\n",
"### Best estimates and confidence intervals\n",
"\n",
"The posterior distribution from our Bayesian data analysis is the key quantity that encodes our inference about the values of the model parameters, given the data and the relevant background information. Often, however, we wish to summarize this result with just a few numbers: the best estimate and a measure of its reliability. \n",
"\n",
"There are a few different options for this. The choice of the most appropriate one depends mainly on the shape of the posterior distribution:\n",
"\n",
"*Symmetric posterior pdfs*: Since the probability (density) associated with any particular value of the parameter is a measure of how much we believe that it lies in the neighbourhood of that point, our best estimate is given by the maximum of the posterior pdf. If we denote the quantity of interest by $X$, with a posterior pdf $P =p(X|D,I)$, then the best estimate of its value $X_0$ is given by the condition $dP/dX|_{X=X_0}=0$. Strictly speaking, we should also check the sign of the second derivative to ensure that $X_0$ represents a maximum.\n",
"\n",
"To obtain a measure of the reliability of this best estimate, we need to look at the width or spread of the posterior pdf about $X_0$. When considering the behaviour of any function in the neighbourhood of a particular point, it is often helpful to carry out a Taylor series expansion; this is simply a standard tool for (locally) approximating a complicated function by a low-order polynomial. The linear term is zero at the maximum and the quadratic term is often the dominating one determining the width of the posterior pdf. Ignoring all the higher-order terms we arrive at the Gaussian approximation"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"_auto8\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"p(X|D,I) \\approx \\frac{1}{\\sigma\\sqrt{2\\pi}} \\exp \\left[ -\\frac{(x-\\mu)^2}{2\\sigma^2} \\right],\n",
"\\label{_auto8} \\tag{11}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"where the mean $\\mu = X_0$ and the variance $\\sigma = \\left( - \\left. \\frac{d^2L}{dX^2} \\right|_{X_0} \\right)^{-1/2}$, where $L$ is the logarithm of the posterior $P$. Our inference about the quantity of interest is conveyed very concisely, therefore, by the statement $X = X_0 \\pm \\sigma$, and"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"$$\n",
"p(X_0-\\sigma < X < X_0+\\sigma | D,I) = \\int_{X_0-\\sigma}^{X_0+\\sigma} p(X|D,I) dX \\approx 0.67.\n",
"$$\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"*Asymmetric posterior pdfs*: While the maximum of the posterior ($X_0$) can still be regarded as giving the best estimate, the true value is now more likely to be on one side of this rather than the other. Alternatively one can compute the mean value, $\\langle X \\rangle = \\int X p(X|D,I) dX$, although this tends to overemphasise very long tails. The best option is probably a compromise that can be employed when having access to a large sample from the posterior (as provided by an MCMC), namely to give the median of this ensamble.\n",
"\n",
"Furthermore, the concept of an error-bar does not seem appropriate in this case, as it implicitly entails the idea of symmetry. A good way of expressing the reliability with which a parameter can be inferred, for an asymmetric posterior pdf, is rather through a *confidence interval*. Since the area under the posterior pdf between $X_1$ and $X_2$ is proportional to how much we believe that $X$ lies in that range, the shortest interval that encloses 67% of the area represents a sensible measure of the uncertainty of the estimate. Obviously we can choose to provide some other degree-of-belief that we think is relevant for the case at hand. Assuming that the posterior pdf has been normalized, to have unit area, we need to find $X_1$ and $X_2$ such that:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"$$\n",
"p(X_1 < X < X_2 | D,I) = \\int_{X_1}^{X_2} p(X|D,I) dX \\approx 0.67, \n",
"$$\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"where the difference $X_2 - X_1$ is as small as possible. The region $X_1 < X < X_2$ is then called the shortest 67% confidence interval. \n",
"\n",
"*Multimodal posterior pdfs*: We can sometimes obtain posteriors which are multimodal; i.e. contains several disconnected regions with large probabilities. There is no difficulty when one of the maxima is very much larger than the others: we can simply ignore the subsidiary solutions, to a good approximation, and concentrate on the global maximum. The problem arises when there are several maxima of comparable magnitude. What do we now mean by a best estimate, and how should we quantify its reliability? The idea of a best estimate and an error-bar, or even a confidence interval, is merely an attempt to summarize the posterior with just two or three numbers; sometimes this just can’t be done, and so these concepts are not valid. For the bimodal case we might be able to characterize the posterior in terms of a few numbers: two best estimates and their associated error-bars, or disjoint confidence intervals. For a general multimodal pdf, the most honest thing we can do is just display the posterior itself.\n",
"\n",
"### Simple Photon Counts: Best estimates and confidence intervals\n",
"\n",
"To compute these numbers for our example, you would run:"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"sampper=np.percentile(samples, [2.5, 16.5, 50, 83.5, 97.5],axis=0).flatten()\n",
"print(\"\"\"\n",
"F_true = {0}\n",
"Based on {1} measurements the posterior point estimates are:\n",
"...F_est = {2:.0f} +/- {3:.0f}\n",
"or using credible intervals:\n",
"...F_est = {4:.0f} (posterior median) \n",
"...F_est in [{5:.0f}, {6:.0f}] (67% credible interval) \n",
"...F_est in [{7:.0f}, {8:.0f}] (95% credible interval) \"\"\"\\\n",
" .format(F_true, N, np.mean(samples), np.std(samples), \\\n",
" sampper[2], sampper[1], sampper[3], sampper[0], sampper[4]))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"`F_true = 1000` \n",
"`Based on 50 measurements the posterior point estimates are:` \n",
"`...F_est = 998 +/- 4` \n",
"`or using credible intervals:` \n",
"`...F_est = 998 (posterior median)` \n",
"`...F_est in [993, 1002] (67% credible interval)` \n",
"`...F_est in [989, 1006] (95% credible interval)` \n",
"\n",
"In this particular example, the posterior pdf is actually a Gaussian (since it is constructed as a product of Gaussians), and the mean and variance from the quadratic approximation will agree exactly with the frequentist approach.\n",
"\n",
"From this final result you might come away with the impression that the Bayesian method is unnecessarily complicated, and in this case it certainly is. Using an MCMC sampler to characterize a one-dimensional normal distribution is a bit like using the Death Star to destroy a beach ball, but we did this here because it demonstrates an approach that can scale to complicated posteriors in many, many dimensions, and can provide nice results in more complicated situations where an analytic likelihood approach is not possible.\n",
"\n",
"Furthermore, as data and models grow in complexity, the two approaches can diverge greatly. \n",
"\n",
"## Bayesian parameter estimation (multiple parameters, covariance)\n",
"* multidimensional posterior pdf:s\n",
"\n",
"* nuisance parameters (e.g. background subtraction?)\n",
"\n",
"* corner plots, covariance, correlations\n",
"\n",
"* best example?\n",
"\n",
"\n",
"\n",
"## Bayesian model selection\n",
"* Bayesian evidence\n",
"\n",
"* Occam's razor\n",
"\n",
"* Best example? How many spectral lines are there?"
]
}
],
"metadata": {},
"nbformat": 4,
"nbformat_minor": 2
"nbformat_minor": 0
}
doc/pub/Bayesian/ipynb/ipynb-Bayesian-src.tar.gz
BIN
View File
Binary file not shown.
doc/pub/Bayesian/pdf/Bayesian-beamer-handouts2x3.pdf
BIN
View File
Binary file not shown.
doc/pub/Bayesian/pdf/Bayesian-beamer.pdf
BIN
View File
Binary file not shown.
doc/pub/Bayesian/pdf/Bayesian-minted.pdf
BIN
View File
Binary file not shown.
doc/src/Bayesian/Bayesian.do.txt
+252 -4
View File
@@ -101,9 +101,10 @@ $$
!et
The importance of this property to data analysis becomes apparent if we replace $A$ and $B$ by hypothesis($H$) and data($D$):
!bt
$$
p(H|D,I) = \frac{p(D|H,I) p(H|I)}{p(D|I)}.
$$
\begin{align}
p(H|D,I) &= \frac{p(D|H,I) p(H|I)}{p(D|I)}.
label{eq:bayes}
\end{align}
!et
The power of Bayes’ theorem lies in the fact that it relates the quantity of interest, the probability that the hypothesis is true given the data, to the term we have a better chance of being able to assign, the probability that we would have observed the measured data if the hypothesis was true.
!eblock
@@ -173,4 +174,251 @@ $$
* The lighthouse problem (two parameters, MCMC)
* Linear fit with outliers (nuisance parameters)
* ...
!eblock
!eblock
!split
===== Example: Is this a fair coin? =====
Let us begin with the analysis of data from a simple coin-tossing experiment.
Given that we had observed 6 heads in 8 flips, would you think it was a fair coin? By fair, we mean that we would be prepared to lay an even 1 : 1 bet on the outcome of a flip being a head or a tail. If we decide that the coin was fair, the question which follows naturally is how sure are we that this was so; if it was not fair, how unfair do we think it was? Furthermore, if we were to continue collecting data for this particular coin, observing the outcomes of additional flips, how would we update our belief on the fairness of the coin?
A sensible way of formulating this problem is to consider a large number of hypotheses about the range in which the bias-weighting of the coin might lie. If we denote the bias-weighting by $H$, then $H = 0$ and $H = 1$ can represent a coin which produces a tail or a head on every flip, respectively. There is a continuum of possibilities for the value of H between these limits, with $H = 0.5$ indicating a fair coin. Our state of knowledge about the fairness, or the degree of unfairness, of the coin is then completely summarized by specifying how much we believe these various propositions to be true.
Let us perform a computer simulation of a coin-tossing experiment. This provides the data that we will be analysing.
@@@CODE src/coinflipping.py from-to: start import modules@end import modules
@@@CODE src/coinflipping.py from-to: start simulation@end simulation
In the light of this data, our inference about the fairness of this coin is summarized by the conditional pdf: $p(H|D,I)$. This is, of course, shorthand for the limiting case of a continuum of propositions for the value of $H$; that is to say, the probability that $H$ lies in an infinitesimally narrow range is given by $p(H|D,I) dH$.
To estimate this posterior pdf, we need to use Bayes’ theorem (ref{eq:bayes}). We will ignore the denominator $p(D|I)$ as it does not involve bias-weighting explicitly, and it will therefore not affect the shape of the desired pdf. At the end we can evaluate the missing constant subsequently from the normalization condition
!bt
\begin{equation}
\int_0^1 p(H|D,I) dH = 1.
label{eq:coin_posterior_norm}
\end{equation}
!et
The prior pdf, $p(H|I)$, represents what we know about the coin given only the information $I$ that we are dealing with a ‘strange coin’. We could keep a very open mind about the nature of the coin; a simple probability assignment which reflects this is a uniform, or flat, prior
!bt
\begin{equation}
p(H|I) = \left\{ \begin{array}{ll}
1 & 0 \le H \le 1, \\
0 & \mathrm{otherwise}.
\end{array} \right.
label{eq:coin_prior_uniform}
\end{equation}
!et
We will get back later to the choice of prior and its effect on the analysis.
This prior state of knowledge, or ignorance, is modified by the data through the likelihood function $p(D|H,I)$. It is a measure of the chance that we would have obtained the data that we actually observed, if the value of the bias-weighting was given (as known). If, in the conditioning information $I$, we assume that the flips of the coin were independent events, so that the outcome of one did not influence that of another, then the probability of obtaining the data `R heads in N tosses' is given by the binomial distribution (we leave a formal definition of this to a statistics textbook)
!bt
\begin{equation}
p(D|H,I) \propto H^R (1-H)^{N-R}.
\end{equation}
!et
It seems reasonable because $H$ is the chance of obtaining a head on any flip, and there were $R$ of them, and $1-H$ is the corresponding probability for a tail, of which there were $N-R$. We note that this binomial distribution also contains a normalization factor, but we will ignore it since it does not depend explicitly on $H$, the quantity of interest. It will be absorbed by the normalization condition (ref{eq:coin_posterior_norm}).
We perform the setup of this Bayesian framework on the computer.
@@@CODE src/coinflipping.py from-to: start bayesian setup@end bayesian setup
The next step is to confront this setup with the simulated data. To get a feel for the result, it is instructive to see how the posterior pdf evolves as we obtain more and more data pertaining to the coin. The results of such an analyses is shown in Fig. ref{fig:coinflipping}.
@@@CODE src/coinflipping.py from-to: start plotting@end plotting
FIGURE:[fig/coinflipping_fig_1.png, width=500 frac=0.95] The evolution of the posterior pdf for the bias-weighting of a coin, as the number of data available increases. The figure on the top left-hand corner of each panel shows the number of data included in the analysis. label{fig:coinflipping}
The panel in the top left-hand corner shows the posterior pdf for $H$ given no data, i.e., it is the same as the prior pdf of Eq. (ref{eq:coin_prior_uniform}). It indicates that we have no more reason to believe that the coin is fair than we have to think that it is double-headed, double-tailed, or of any other intermediate bias-weighting.
The first flip is obviously tails. At this point we have no evidence that the coin has a side with heads, as indicated by the pdf going to zero as $H \to 1$. The second flip is obviously heads and we have now excluded both extreme options $H=0$ (double-tailed) and $H=1$ (double-headed). We can note that the posterior at this point has the simple form $p(H|D,I) = H(1-H)$ for $0 \le H \le 1$.
The remainder of Fig. ref{fig:coinflipping} shows how the posterior pdf evolves as the number of data analysed becomes larger and larger. We see that the position of the maximum moves around, but that the amount by which it does so decreases with the increasing number of observations. The width of the posterior pdf also becomes narrower with more data, indicating that we are becoming increasingly confident in our estimate of the bias-weighting. For the coin in this example, the best estimate of $H$ eventually converges to 0.6, which, of course, was the value chosen to simulate the flips.
!split
===== A few words on different priors =====
* uniform
* Gaussian
* Jeffrey's prior
Repeat the coin flipping experiment with other priors.
!split
===== Bayesian parameter estimation (single parameter) =====
!bblock
We will now consider the very important task of model parameter estimation using statistical inference.
[CF: maybe stress that model parameters are not random variables, and the meaning of parameter estimation is therefore very different between frequentist and bayesian approaches.]
Throughout this section we will consider a specific example that involves a model with a single parameter: ``Measured flux from a star''.
!eblock
!split
=== Example: Measured flux from a star ===
Adapted from the blog "Pythonic Perambulations": "http://jakevdp.github.io" by Jake VanderPlas.
Imagine that we point our telescope to the sky, and observe the light coming from a single star. For the time being, we'll assume that the star's true flux is constant with time, i.e. that is it has a fixed value $F_\mathrm{true}$ (we'll also ignore effects like sky noise and other sources of systematic error). We'll assume that we perform a series of $N$ measurements with our telescope, where the ith measurement reports the observed photon flux $F_i$ and error $e_i$[^errors].
The question is, given this set of measurements $D = \{F_i, e_i\}$, what is our best estimate of the true flux $F_\mathrm{true}$?
[^errors]: We'll make the reasonable assumption that errors are Gaussian. In a Frequentist perspective, $e_i$ is the standard deviation of the results of a single measurement event in the limit of repetitions of *that event*. In the Bayesian perspective, $e_i$ is the standard deviation of the (Gaussian) probability distribution describing our knowledge of that particular measurement given its observed value.
Because the measurements are number counts, a Poisson distribution is a good approximation to the measurement process:
@@@CODE src/singlephotoncount.py from-to: start generate data@end generate data
Now let's make a simple visualization of the ``observed'' data, see Fig. ref{fig:flux}.
@@@CODE src/singlephotoncount.py from-to: start visualize data@end visualize data
FIGURE:[fig/singlephotoncount_fig_1.png, width=400 frac=0.8] Single photon counts (flux measurements). label{fig:flux}
These measurements each have a different error $e_i$ which is estimated from Poisson statistics using the standard square-root rule. In this toy example we already know the true flux $F_\mathrm{true}$, but the question is this: given our measurements and errors, what is our best estimate of the true flux?
Let's take a look at the frequentist and Bayesian approaches to solving this.
=== Simple Photon Counts: Frequentist Approach ===
We'll start with the classical frequentist maximum likelihood approach. Given a single observation $D_i = (F_i, e_i)$, we can compute the probability distribution of the measurement given the true flux Ftrue given our assumption of Gaussian errors
!bt
\begin{equation}
p(D_i | F_\mathrm{true}, I) = \frac{1}{\sqrt{2\pi e_i^2}} \exp \left( \frac{-(F_i-F_\mathrm{true})^2}{2e_i^2} \right).
\end{equation}
!et
This should be read ``the probability of $D_i$ given $F_\mathrm{true}$
equals ...''. You should recognize this as a normal distribution with mean $F_\mathrm{true}$ and standard deviation $e_i$.
We construct the *likelihood function* by computing the product of the probabilities for each data point
!bt
\begin{equation}
\mathcal{L}(D | F_\mathrm{true}, I) = \prod_{i=1}^N p(D_i | F_\mathrm{true}, I),
\end{equation}
!et
here $D = \{D_i\}$ represents the entire set of measurements. Because the value of the likelihood can become very small, it is often more convenient to instead compute the log-likelihood. Combining the previous two equations and computing the log, we have
!bt
\begin{equation}
\log\mathcal{L} = -\frac{1}{2} \sum_{i=1}^N \left[ \log(2\pi e_i^2) + \frac{(F_i-F_\mathrm{true})^2}{e_i^2} \right].
\end{equation}
!et
What we'd like to do is determine $F_\mathrm{true}$ such that the likelihood is maximized. For this simple problem, the maximization can be computed analytically (i.e. by setting $d\log\mathcal{L}/d F_\mathrm{true} = 0$). This results in the following observed estimate of $F_\mathrm{true}$
!bt
\begin{equation}
F_\mathrm{est} = \sum_{i=1}^N w_i F_i; \quad w_i = 1/e_i^2.
\end{equation}
!et
Notice that in the special case of all errors $e_i$ being equal, this reduces to
!bt
\begin{equation}
F_\mathrm{est} = \frac{1}{N} \sum_{i=1} F_i.
\end{equation}
!et
That is, in agreement with intuition, $F_\mathrm{est}$ is simply the mean of the observed data when errors are equal.
We can go further and ask what the error of our estimate is. In the frequentist approach, this can be accomplished by fitting a Gaussian approximation to the likelihood curve at maximum; in this simple case this can also be solved analytically (the sum of Gaussians is also a Gaussian). It can be shown that the standard deviation of this Gaussian approximation is
!bt
\begin{equation}
\sigma_\mathrm{est} = \sum_{i=1}^N w_i.
\end{equation}
!et
These results are fairly simple calculations; let's evaluate them for our toy dataset:
@@@CODE src/singlephotoncount.py from-to: start frequentist@end frequentist
`F_true = 1000` <linebreak>
`F_est = 998 +/- 4 (based on 50 measurements)` <linebreak>
We find that for 50 measurements of the flux, our estimate has an error of about 0.4% and is consistent with the input value.
=== Simple Photon Counts: Bayesian Approach ===
The Bayesian approach, as you might expect, begins and ends with probabilities. Our hypothesis is that the star has a constant flux $F_\mathrm{true}$. It recognizes that what we fundamentally want to compute is our knowledge of the parameters in question given the data and other information (such as our knowledge of uncertainties for the observed values), i.e. in this case, $p(F_\mathrm{true} | D,I)$.
Note that this formulation of the problem is fundamentally contrary to the frequentist philosophy, which says that probabilities have no meaning for model parameters like $F_\mathrm{true}$. Nevertheless, within the Bayesian philosophy this is perfectly acceptable.
To compute this result, Bayesians next apply Bayes' Theorem (ref{eq:bayes}).
If we set the prior $p(F_\mathrm{true}|I) \propto 1$ (a flat prior), we find
$p(F_\mathrm{true}|D,I) \propto p(D | F_\mathrm{true},I) \equiv \mathcal{L}(D | F_\mathrm{true},I)$
and the Bayesian probability is maximized at precisely the same value as the frequentist result! So despite the philosophical differences, we see that (for this simple problem at least) the Bayesian and frequentist point estimates are equivalent.
=== A note about priors ===
The prior allows inclusion of other information into the computation, which becomes very useful in cases where multiple measurement strategies are being combined to constrain a single model. The necessity to specify a prior, however, is one of the more controversial pieces of Bayesian analysis.
A frequentist will point out that the prior is problematic when no true prior information is available. Though it might seem straightforward to use a noninformative prior like the flat prior mentioned above, there are some "surprisingly subtleties": "http://normaldeviate.wordpress.com/2013/07/13/lost-causes-in-statistics-ii-noninformative- priors/comment-page-1/" involved. It turns out that in many situations, a truly noninformative prior does not exist! Frequentists point out that the subjective choice of a prior which necessarily biases your result has no place in statistical data analysis.
A Bayesian would counter that frequentism doesn't solve this problem, but simply skirts the question. Frequentism can often be viewed as simply a special case of the Bayesian approach for some (implicit) choice of the prior: a Bayesian would say that it's better to make this implicit choice explicit, even if the choice might include some subjectivity.
=== Simple Photon Counts: Bayesian approach in practice ===
Leaving these philosophical debates aside for the time being, let's address how Bayesian results are generally computed in practice. For a one parameter problem like the one considered here, it's as simple as computing the posterior probability $p(F_\mathrm{true} | D,I)$ as a function of $F_\mathrm{true}$: this is the distribution reflecting our knowledge of the parameter $F_\mathrm{true}$.
But as the dimension of the model grows, this direct approach becomes increasingly intractable. For this reason, Bayesian calculations often depend on sampling methods such as Markov Chain Monte Carlo (MCMC). For this practical example, let us apply an MCMC approach using Dan Foreman-Mackey's "emcee": "http://dan.iel.fm/emcee/current/" package. Keep in mind here that the goal is to generate a set of points drawn from the posterior probability distribution, and to use those points to determine the answer we seek.
To perform this MCMC, we start by defining Python functions for the prior $p(F_\mathrm{true} | I)$, the likelihood $p(D | F_\mathrm{true},I)$, and the posterior $p(F_\mathrm{true} | D,I)$, noting that none of these need be properly normalized. Our model here is one-dimensional, but to handle multi-dimensional models we'll define the model in terms of an array of parameters $\boldsymbol{\alpha}$, which in this case is $\boldsymbol{\alpha} = [F_\mathrm{true}]$
@@@CODE src/singlephotoncount.py from-to: start bayesian setup@end bayesian setup
Now we set up the problem, including generating some random starting guesses for the multiple chains of points.
@@@CODE src/singlephotoncount.py from-to: start bayesian mcmc@end bayesian mcmc
If this all worked correctly, the array sample should contain a series of 50,000 points drawn from the posterior. Let's plot them and check. See results in Fig. ref{fig:flux-bayesian}.
@@@CODE src/singlephotoncount.py from-to: start visualize bayesian@end visualize bayesian
FIGURE:[fig/singlephotoncount_fig_2.png, width=400 frac=0.8] Bayesian posterior pdf (represented by a histogram of MCMC samples) from flux measurements. label{fig:flux-bayesian}
=== Best estimates and confidence intervals ===
The posterior distribution from our Bayesian data analysis is the key quantity that encodes our inference about the values of the model parameters, given the data and the relevant background information. Often, however, we wish to summarize this result with just a few numbers: the best estimate and a measure of its reliability.
There are a few different options for this. The choice of the most appropriate one depends mainly on the shape of the posterior distribution:
*Symmetric posterior pdfs*: Since the probability (density) associated with any particular value of the parameter is a measure of how much we believe that it lies in the neighbourhood of that point, our best estimate is given by the maximum of the posterior pdf. If we denote the quantity of interest by $X$, with a posterior pdf $P =p(X|D,I)$, then the best estimate of its value $X_0$ is given by the condition $dP/dX|_{X=X_0}=0$. Strictly speaking, we should also check the sign of the second derivative to ensure that $X_0$ represents a maximum.
To obtain a measure of the reliability of this best estimate, we need to look at the width or spread of the posterior pdf about $X_0$. When considering the behaviour of any function in the neighbourhood of a particular point, it is often helpful to carry out a Taylor series expansion; this is simply a standard tool for (locally) approximating a complicated function by a low-order polynomial. The linear term is zero at the maximum and the quadratic term is often the dominating one determining the width of the posterior pdf. Ignoring all the higher-order terms we arrive at the Gaussian approximation
!bt
\begin{equation}
p(X|D,I) \approx \frac{1}{\sigma\sqrt{2\pi}} \exp \left[ -\frac{(x-\mu)^2}{2\sigma^2} \right],
\end{equation}
!et
where the mean $\mu = X_0$ and the variance $\sigma = \left( - \left. \frac{d^2L}{dX^2} \right|_{X_0} \right)^{-1/2}$, where $L$ is the logarithm of the posterior $P$. Our inference about the quantity of interest is conveyed very concisely, therefore, by the statement $X = X_0 \pm \sigma$, and
!bt
$$
p(X_0-\sigma < X < X_0+\sigma | D,I) = \int_{X_0-\sigma}^{X_0+\sigma} p(X|D,I) dX \approx 0.67.
$$
!et
*Asymmetric posterior pdfs*: While the maximum of the posterior ($X_0$) can still be regarded as giving the best estimate, the true value is now more likely to be on one side of this rather than the other. Alternatively one can compute the mean value, $\langle X \rangle = \int X p(X|D,I) dX$, although this tends to overemphasise very long tails. The best option is probably a compromise that can be employed when having access to a large sample from the posterior (as provided by an MCMC), namely to give the median of this ensamble.
Furthermore, the concept of an error-bar does not seem appropriate in this case, as it implicitly entails the idea of symmetry. A good way of expressing the reliability with which a parameter can be inferred, for an asymmetric posterior pdf, is rather through a *confidence interval*. Since the area under the posterior pdf between $X_1$ and $X_2$ is proportional to how much we believe that $X$ lies in that range, the shortest interval that encloses 67% of the area represents a sensible measure of the uncertainty of the estimate. Obviously we can choose to provide some other degree-of-belief that we think is relevant for the case at hand. Assuming that the posterior pdf has been normalized, to have unit area, we need to find $X_1$ and $X_2$ such that:
!bt
$$
p(X_1 < X < X_2 | D,I) = \int_{X_1}^{X_2} p(X|D,I) dX \approx 0.67,
$$
!et
where the difference $X_2 - X_1$ is as small as possible. The region $X_1 < X < X_2$ is then called the shortest 67% confidence interval.
*Multimodal posterior pdfs*: We can sometimes obtain posteriors which are multimodal; i.e. contains several disconnected regions with large probabilities. There is no difficulty when one of the maxima is very much larger than the others: we can simply ignore the subsidiary solutions, to a good approximation, and concentrate on the global maximum. The problem arises when there are several maxima of comparable magnitude. What do we now mean by a best estimate, and how should we quantify its reliability? The idea of a best estimate and an error-bar, or even a confidence interval, is merely an attempt to summarize the posterior with just two or three numbers; sometimes this just can’t be done, and so these concepts are not valid. For the bimodal case we might be able to characterize the posterior in terms of a few numbers: two best estimates and their associated error-bars, or disjoint confidence intervals. For a general multimodal pdf, the most honest thing we can do is just display the posterior itself.
=== Simple Photon Counts: Best estimates and confidence intervals ===
To compute these numbers for our example, you would run:
@@@CODE src/singlephotoncount.py from-to: start bayesian CI@end bayesian CI
`F_true = 1000` <linebreak>
`Based on 50 measurements the posterior point estimates are:` <linebreak>
`...F_est = 998 +/- 4` <linebreak>
`or using credible intervals:` <linebreak>
`...F_est = 998 (posterior median)` <linebreak>
`...F_est in [993, 1002] (67% credible interval)` <linebreak>
`...F_est in [989, 1006] (95% credible interval)` <linebreak>
In this particular example, the posterior pdf is actually a Gaussian (since it is constructed as a product of Gaussians), and the mean and variance from the quadratic approximation will agree exactly with the frequentist approach.
From this final result you might come away with the impression that the Bayesian method is unnecessarily complicated, and in this case it certainly is. Using an MCMC sampler to characterize a one-dimensional normal distribution is a bit like using the Death Star to destroy a beach ball, but we did this here because it demonstrates an approach that can scale to complicated posteriors in many, many dimensions, and can provide nice results in more complicated situations where an analytic likelihood approach is not possible.
Furthermore, as data and models grow in complexity, the two approaches can diverge greatly.
!split
===== Bayesian parameter estimation (multiple parameters, covariance) =====
!bblock
* multidimensional posterior pdf:s
* nuisance parameters (e.g. background subtraction?)
* corner plots, covariance, correlations
* best example?
!eblock
!split
===== Bayesian model selection =====
!bblock
* Bayesian evidence
* Occam's razor
* Best example? How many spectral lines are there?
!eblock
doc/src/Bayesian/fig/coinflipping_fig_1.png
BIN
View File
Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 42 KiB

doc/src/Bayesian/fig/singlephotoncount_fig_1.png
BIN
View File
Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 19 KiB

doc/src/Bayesian/fig/singlephotoncount_fig_2.png
BIN
View File
Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 12 KiB

doc/src/Bayesian/fig/sivia_fig_1_2.png
BIN
View File
Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 18 KiB

doc/src/Bayesian/makepdf.sh
Executable
+48
View File
@@ -0,0 +1,48 @@
#!/bin/sh
set -x
function system {
"$@"
if [ $? -ne 0 ]; then
echo "make.sh: unsuccessful command $@"
echo "abort!"
exit 1
fi
}
if [ $# -eq 0 ]; then
echo 'bash make.sh slides1|slides2'
exit 1
fi
name=$1
rm -f *.tar.gz
opt="--encoding=utf-8"
# Note: Makefile examples contain constructions like ${PROG} which
# looks like Mako constructions, but they are not. Use --no_mako
# to turn off Mako processing.
opt="--no_mako"
rm -f *.aux
# Ordinary plain LaTeX document
rm -f *.aux # important after beamer
system doconce format pdflatex $name --minted_latex_style=trac --latex_admon=paragraph $opt
system doconce ptex2tex $name envir=minted
# Add special packages
doconce subst "% Add user's preamble" "\g<1>\n\\usepackage{simplewick}" $name.tex
doconce replace 'section{' 'section*{' $name.tex
pdflatex -shell-escape $name
pdflatex -shell-escape $name
mv -f $name.pdf ${name}-minted.pdf
cp $name.tex ${name}-plain-minted.tex
# Publish
dest=../../pub
if [ ! -d $dest/$name ]; then
mkdir $dest/$name
mkdir $dest/$name/pdf
fi
cp ${name}*.pdf $dest/$name/pdf
doc/src/Bayesian/src/coinflipping.py
+49
View File
@@ -0,0 +1,49 @@
# start import modules
import numpy as np
import matplotlib.pyplot as plt
# end import modules
savefig=True
# start simulation
np.random.seed(999) # for reproducibility
a=0.6 # biased coin
flips=np.random.rand(2**12) # simulates 4096 coin flips
heads=flips<a # boolean array, heads[i]=True if flip i is heads
# end simulation
# start bayesian setup
def prior(H):
p=np.zeros_like(H)
p[(0<=x)&(x<=1)]=1 # allowed range: 0<=H<=1
return p # uniform prior
def likelihood(H,data):
N = len(data)
no_of_heads = sum(data)
no_of_tails = N - no_of_heads
return H**no_of_heads * (1-H)**no_of_tails
def posterior(H,data):
p=prior(H)*likelihood(H,data)
norm=np.trapz(p,H)
return p/norm
# end bayesian setup
# start plotting
x=np.linspace(0,1,100)
fig, axs = plt.subplots(nrows=4,ncols=3,sharex=True,sharey='row')
axs_vec=np.reshape(axs,-1)
axs_vec[0].plot(x,prior(x))
for ndouble in range(11):
ax=axs_vec[1+ndouble]
ax.plot(x,posterior(x,heads[:2**ndouble]))
ax.text(0.1, 0.8, '$N={0}$'.format(2**ndouble), transform=ax.transAxes)
for row in range(4): axs[row,0].set_ylabel('$p(H|D_\mathrm{obs},I)$')
for col in range(3): axs[-1,col].set_xlabel('$H$')
# end plotting
if savefig:
fig.savefig('../fig/coinflipping_fig_1.png')
if not savefig:
plt.show()
doc/src/Bayesian/src/environment.yml
+10
View File
@@ -0,0 +1,10 @@
name: ML-bayesian-env
channels:
- conda-forge
dependencies:
- python=3.6
- numpy
- matplotlib
- emcee
- gpy # Gaussian Processes
- corner # Plot tool (for e.g. projected posterior distributions)
doc/src/Bayesian/src/singlephotoncount.py
+94
View File
@@ -0,0 +1,94 @@
# start import modules
import numpy as np
from scipy import stats
import matplotlib.pyplot as plt
import emcee
# end import modules
savefig=True
# start generate data
np.random.seed(1) # for repeatability
F_true = 1000 # true flux, say number of photons measured in 1 second
N = 50 # number of measurements
F = stats.poisson(F_true).rvs(N)
# N measurements of the flux
e = np.sqrt(F) # errors on Poisson counts estimated via square root
# end generate data
# start visualize data
fig, ax = plt.subplots()
ax.errorbar(F, np.arange(N), xerr=e, fmt='ok', ecolor='gray', alpha=0.5)
ax.vlines([F_true], 0, N, linewidth=5, alpha=0.2)
ax.set_xlabel("Flux");ax.set_ylabel("measurement number");
# end visualize data
if savefig:
fig.savefig('../fig/singlephotoncount_fig_1.png')
# start frequentist
w=1./e**2
print("""
F_true = {0}
F_est = {1:.0f} +/- {2:.0f} (based on {3} measurements) """\
.format(F_true, (w * F).sum() / w.sum(), w.sum() ** -0.5, N))
# end frequentist
# start bayesian setup
def log_prior(alpha):
return 0 # flat prior
def log_likelihood(alpha, F, e):
return -0.5 * np.sum(np.log(2 * np.pi * e ** 2) \
+ (F - alpha[0]) ** 2 / e ** 2)
def log_posterior(alpha, F, e):
return log_prior(alpha) + log_likelihood(alpha, F, e)
# end bayesian setup
# start bayesian mcmc
ndim = 1 # number of parameters in the model
nwalkers = 50 # number of MCMC walkers
nburn = 1000 # "burn-in" period to let chains stabilize
nsteps = 2000 # number of MCMC steps to take
# we'll start at random locations between 0 and 2000
starting_guesses = 2000 * np.random.rand(nwalkers, ndim)
sampler = emcee.EnsembleSampler(nwalkers, ndim, log_posterior, args=[F,e])
sampler.run_mcmc(starting_guesses, nsteps)
# Shape of sampler.chain = (nwalkers, nsteps, ndim)
# Flatten the sampler chain and discard burn-in points:
samples = sampler.chain[:, nburn:, :].reshape((-1, ndim))
# end bayesian mcmc
# start visualize bayesian
fig, ax = plt.subplots()
ax.hist(samples, bins=50, histtype="stepfilled", alpha=0.3, normed=True)
ax.set_xlabel(r'$F_\mathrm{est}$')
ax.set_ylabel(r'$p(F_\mathrm{est}|D,I)$')
# end visualize bayesian
if savefig:
fig.savefig('../fig/singlephotoncount_fig_2.png')
# plot a best-fit Gaussian
F_est = np.linspace(975, 1025)
pdf = stats.norm(np.mean(samples), np.std(samples)).pdf(F_est)
ax.plot(F_est, pdf, '-k')
# start bayesian CI
sampper=np.percentile(samples, [2.5, 16.5, 50, 83.5, 97.5],axis=0).flatten()
print("""
F_true = {0}
Based on {1} measurements the posterior point estimates are:
...F_est = {2:.0f} +/- {3:.0f}
or using credible intervals:
...F_est = {4:.0f} (posterior median)
...F_est in [{5:.0f}, {6:.0f}] (67% credible interval)
...F_est in [{7:.0f}, {8:.0f}] (95% credible interval) """\
.format(F_true, N, np.mean(samples), np.std(samples), \
sampper[2], sampper[1], sampper[3], sampper[0], sampper[4]))
# end bayesian CI
if not savefig:
plt.show()