
<!DOCTYPE html>

<html lang="en">
  <head>
    <meta charset="utf-8" />
    <meta name="viewport" content="width=device-width, initial-scale=1.0" /><meta name="generator" content="Docutils 0.18.1: http://docutils.sourceforge.net/" />

    <title>Case Study: Outliers and Bayesian Integrals &#8212; lsqfit 13.0.1 documentation</title>
    <link rel="stylesheet" type="text/css" href="_static/pygments.css" />
    <link rel="stylesheet" type="text/css" href="_static/pyramid.css" />
    <script data-url_root="./" id="documentation_options" src="_static/documentation_options.js"></script>
    <script src="_static/jquery.js"></script>
    <script src="_static/underscore.js"></script>
    <script src="_static/_sphinx_javascript_frameworks_compat.js"></script>
    <script src="_static/doctools.js"></script>
    <link rel="index" title="Index" href="genindex.html" />
    <link rel="search" title="Search" href="search.html" />
    <link rel="next" title="lsqfit - Nonlinear Least Squares Fitting" href="lsqfit.html" />
    <link rel="prev" title="Case Study: Fitting a Spline" href="case-spline.html" />
<link rel="stylesheet" href="https://fonts.googleapis.com/css?family=Neuton&amp;subset=latin" type="text/css" media="screen" charset="utf-8" />
<link rel="stylesheet" href="https://fonts.googleapis.com/css?family=Nobile:regular,italic,bold,bolditalic&amp;subset=latin" type="text/css" media="screen" charset="utf-8" />
<!--[if lte IE 6]>
<link rel="stylesheet" href="_static/ie6.css" type="text/css" media="screen" charset="utf-8" />
<![endif]-->

  </head><body>

    <div class="related" role="navigation" aria-label="related navigation">
      <h3>Navigation</h3>
      <ul>
        <li class="right" style="margin-right: 10px">
          <a href="genindex.html" title="General Index"
             accesskey="I">index</a></li>
        <li class="right" >
          <a href="py-modindex.html" title="Python Module Index"
             >modules</a> |</li>
        <li class="right" >
          <a href="lsqfit.html" title="lsqfit - Nonlinear Least Squares Fitting"
             accesskey="N">next</a> |</li>
        <li class="right" >
          <a href="case-spline.html" title="Case Study: Fitting a Spline"
             accesskey="P">previous</a> |</li>
        <li class="nav-item nav-item-0"><a href="index.html">lsqfit 13.0.1 documentation</a> &#187;</li>
        <li class="nav-item nav-item-this"><a href="">Case Study: Outliers and Bayesian Integrals</a></li> 
      </ul>
    </div>  

    <div class="document">
      <div class="documentwrapper">
        <div class="bodywrapper">
          <div class="body" role="main">
            
  <section id="case-study-outliers-and-bayesian-integrals">
<span id="outliers"></span><h1>Case Study: Outliers and Bayesian Integrals<a class="headerlink" href="#case-study-outliers-and-bayesian-integrals" title="Permalink to this heading">¶</a></h1>
<p>In this case study, we analyze a fit with outliers in the data that
distort the least-squares solution. We show one approach to dealing with
the outliers that requires using Bayesian integrals
in place of least-squares fitting,
to fit the data while also modeling the outliers.</p>
<p>This case study is adapted from an example by Jake Vanderplas
on his <a class="reference external" href="http://jakevdp.github.io/blog/2014/06/06/frequentism-and-bayesianism-2-when-results-differ">Python blog</a>.
It is also discussed in the documentation for the <code class="xref py py-mod docutils literal notranslate"><span class="pre">vegas</span></code> module.</p>
<section id="the-problem">
<h2>The Problem<a class="headerlink" href="#the-problem" title="Permalink to this heading">¶</a></h2>
<p>We want to extrapolate a set of data values <code class="docutils literal notranslate"><span class="pre">y</span></code> to <code class="docutils literal notranslate"><span class="pre">x=0</span></code> fitting
a linear fit function (<code class="docutils literal notranslate"><span class="pre">fitfcn(x,p)</span></code>) to the data:</p>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span><span class="kn">import</span> <span class="nn">numpy</span> <span class="k">as</span> <span class="nn">np</span>

<span class="kn">import</span> <span class="nn">gvar</span> <span class="k">as</span> <span class="nn">gv</span>
<span class="kn">import</span> <span class="nn">lsqfit</span>

<span class="k">def</span> <span class="nf">main</span><span class="p">():</span>
    <span class="c1"># least-squares fit to the data</span>
    <span class="n">x</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">array</span><span class="p">([</span>
        <span class="mf">0.2</span><span class="p">,</span> <span class="mf">0.4</span><span class="p">,</span> <span class="mf">0.6</span><span class="p">,</span> <span class="mf">0.8</span><span class="p">,</span> <span class="mf">1.</span><span class="p">,</span>
        <span class="mf">1.2</span><span class="p">,</span> <span class="mf">1.4</span><span class="p">,</span> <span class="mf">1.6</span><span class="p">,</span> <span class="mf">1.8</span><span class="p">,</span> <span class="mf">2.</span><span class="p">,</span>
        <span class="mf">2.2</span><span class="p">,</span> <span class="mf">2.4</span><span class="p">,</span> <span class="mf">2.6</span><span class="p">,</span> <span class="mf">2.8</span><span class="p">,</span> <span class="mf">3.</span><span class="p">,</span>
        <span class="mf">3.2</span><span class="p">,</span> <span class="mf">3.4</span><span class="p">,</span> <span class="mf">3.6</span><span class="p">,</span> <span class="mf">3.8</span>
        <span class="p">])</span>
    <span class="n">y</span> <span class="o">=</span> <span class="n">gv</span><span class="o">.</span><span class="n">gvar</span><span class="p">([</span>
        <span class="s1">&#39;0.38(20)&#39;</span><span class="p">,</span> <span class="s1">&#39;2.89(20)&#39;</span><span class="p">,</span> <span class="s1">&#39;0.85(20)&#39;</span><span class="p">,</span> <span class="s1">&#39;0.59(20)&#39;</span><span class="p">,</span> <span class="s1">&#39;2.88(20)&#39;</span><span class="p">,</span>
        <span class="s1">&#39;1.44(20)&#39;</span><span class="p">,</span> <span class="s1">&#39;0.73(20)&#39;</span><span class="p">,</span> <span class="s1">&#39;1.23(20)&#39;</span><span class="p">,</span> <span class="s1">&#39;1.68(20)&#39;</span><span class="p">,</span> <span class="s1">&#39;1.36(20)&#39;</span><span class="p">,</span>
        <span class="s1">&#39;1.51(20)&#39;</span><span class="p">,</span> <span class="s1">&#39;1.73(20)&#39;</span><span class="p">,</span> <span class="s1">&#39;2.16(20)&#39;</span><span class="p">,</span> <span class="s1">&#39;1.85(20)&#39;</span><span class="p">,</span> <span class="s1">&#39;2.00(20)&#39;</span><span class="p">,</span>
        <span class="s1">&#39;2.11(20)&#39;</span><span class="p">,</span> <span class="s1">&#39;2.75(20)&#39;</span><span class="p">,</span> <span class="s1">&#39;0.86(20)&#39;</span><span class="p">,</span> <span class="s1">&#39;2.73(20)&#39;</span>
        <span class="p">])</span>
    <span class="n">fit</span> <span class="o">=</span> <span class="n">lsqfit</span><span class="o">.</span><span class="n">nonlinear_fit</span><span class="p">(</span><span class="n">data</span><span class="o">=</span><span class="p">(</span><span class="n">x</span><span class="p">,</span> <span class="n">y</span><span class="p">),</span> <span class="n">prior</span><span class="o">=</span><span class="n">make_prior</span><span class="p">(),</span> <span class="n">fcn</span><span class="o">=</span><span class="n">fitfcn</span><span class="p">)</span>
    <span class="nb">print</span><span class="p">(</span><span class="n">fit</span><span class="p">)</span>

<span class="k">def</span> <span class="nf">fitfcn</span><span class="p">(</span><span class="n">x</span><span class="p">,</span> <span class="n">p</span><span class="p">):</span>
    <span class="n">c</span> <span class="o">=</span> <span class="n">p</span><span class="p">[</span><span class="s1">&#39;c&#39;</span><span class="p">]</span>
    <span class="k">return</span> <span class="n">c</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span> <span class="o">+</span> <span class="n">c</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span> <span class="o">*</span> <span class="n">x</span>

<span class="k">def</span> <span class="nf">make_prior</span><span class="p">():</span>
    <span class="n">prior</span> <span class="o">=</span> <span class="n">gv</span><span class="o">.</span><span class="n">BufferDict</span><span class="p">(</span><span class="n">c</span><span class="o">=</span><span class="n">gv</span><span class="o">.</span><span class="n">gvar</span><span class="p">([</span><span class="s1">&#39;0(5)&#39;</span><span class="p">,</span> <span class="s1">&#39;0(5)&#39;</span><span class="p">]))</span>
    <span class="k">return</span> <span class="n">prior</span>


<span class="k">if</span> <span class="vm">__name__</span> <span class="o">==</span> <span class="s1">&#39;__main__&#39;</span><span class="p">:</span>
    <span class="n">main</span><span class="p">()</span>
</pre></div>
</div>
<p>The fit is not good, with a <code class="docutils literal notranslate"><span class="pre">chi**2</span></code> per degree of freedom that is
much larger than one, despite rather broad priors for the intercept and
slope:</p>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span><span class="n">Least</span> <span class="n">Square</span> <span class="n">Fit</span><span class="p">:</span>
  <span class="n">chi2</span><span class="o">/</span><span class="n">dof</span> <span class="p">[</span><span class="n">dof</span><span class="p">]</span> <span class="o">=</span> <span class="mi">13</span> <span class="p">[</span><span class="mi">19</span><span class="p">]</span>    <span class="n">Q</span> <span class="o">=</span> <span class="mf">1.2e-40</span>    <span class="n">logGBF</span> <span class="o">=</span> <span class="o">-</span><span class="mf">117.45</span>

<span class="n">Parameters</span><span class="p">:</span>
            <span class="n">c</span> <span class="mi">0</span>   <span class="mf">1.149</span> <span class="p">(</span><span class="mi">95</span><span class="p">)</span>     <span class="p">[</span>  <span class="mf">0.0</span> <span class="p">(</span><span class="mf">5.0</span><span class="p">)</span> <span class="p">]</span>  
              <span class="mi">1</span>   <span class="mf">0.261</span> <span class="p">(</span><span class="mi">42</span><span class="p">)</span>     <span class="p">[</span>  <span class="mf">0.0</span> <span class="p">(</span><span class="mf">5.0</span><span class="p">)</span> <span class="p">]</span>  

<span class="n">Settings</span><span class="p">:</span>
  <span class="n">svdcut</span><span class="o">/</span><span class="n">n</span> <span class="o">=</span> <span class="mf">1e-12</span><span class="o">/</span><span class="mi">0</span>    <span class="n">tol</span> <span class="o">=</span> <span class="p">(</span><span class="mf">1e-08</span><span class="o">*</span><span class="p">,</span><span class="mf">1e-10</span><span class="p">,</span><span class="mf">1e-10</span><span class="p">)</span>    <span class="p">(</span><span class="n">itns</span><span class="o">/</span><span class="n">time</span> <span class="o">=</span> <span class="mi">5</span><span class="o">/</span><span class="mf">0.0</span><span class="p">)</span>

</pre></div>
</div>
<p>The problem is evident if we plot the data:</p>
<a class="reference internal image-reference" href="_images/case-outliers1.png"><img alt="_images/case-outliers1.png" src="_images/case-outliers1.png" style="width: 60%;" /></a>
<p>At least three of the data points are outliers: they disagree with other
nearby points by several standard deviations. These outliers have a big
impact on the fit (dashed line). In particular they pull the <code class="docutils literal notranslate"><span class="pre">x=0</span></code> intercept (<code class="docutils literal notranslate"><span class="pre">fit.p['c'][0]</span></code>)
up above one, while the rest of the data suggest an intercept of 0.5
or less.</p>
</section>
<section id="a-solution">
<h2>A Solution<a class="headerlink" href="#a-solution" title="Permalink to this heading">¶</a></h2>
<p>There are many <em>ad hoc</em> prescriptions for handling outliers. In the best
of situations one would have an explanation for the outliers and seek
to model them accordingly. For example,
we might know that some fraction <code class="docutils literal notranslate"><span class="pre">w</span></code> of the time our detector
malfunctions, resulting in much larger measurement errors than usual.
This model can be represented by a more complicated probability
density function (PDF) for the data that consists of a linear combination
of the normal PDF with another PDF that is similar but with much larger
errors. The relative weights assigned to these two terms would be <code class="docutils literal notranslate"><span class="pre">1-w</span></code>
and <code class="docutils literal notranslate"><span class="pre">w</span></code>, respectively.</p>
<p>A modified data prior of this sort is incompatible with the least-squares
code in <a class="reference internal" href="lsqfit.html#module-lsqfit" title="lsqfit: Nonlinear least squares fitting."><code class="xref py py-mod docutils literal notranslate"><span class="pre">lsqfit</span></code></a>. Here we will incorporate it by replacing the
least-squares analysis with a Bayesian integral, where the normal PDF is
replaced by a modified PDF of the sort described above.
The complete code for this analysis is
as follows:</p>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span><span class="kn">import</span> <span class="nn">numpy</span> <span class="k">as</span> <span class="nn">np</span>

<span class="kn">import</span> <span class="nn">gvar</span> <span class="k">as</span> <span class="nn">gv</span>
<span class="kn">import</span> <span class="nn">lsqfit</span>
<span class="kn">import</span> <span class="nn">vegas</span>

<span class="k">def</span> <span class="nf">main</span><span class="p">():</span>
    <span class="c1">### 1) least-squares fit to the data</span>
    <span class="n">x</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">array</span><span class="p">([</span>
        <span class="mf">0.2</span><span class="p">,</span> <span class="mf">0.4</span><span class="p">,</span> <span class="mf">0.6</span><span class="p">,</span> <span class="mf">0.8</span><span class="p">,</span> <span class="mf">1.</span><span class="p">,</span>
        <span class="mf">1.2</span><span class="p">,</span> <span class="mf">1.4</span><span class="p">,</span> <span class="mf">1.6</span><span class="p">,</span> <span class="mf">1.8</span><span class="p">,</span> <span class="mf">2.</span><span class="p">,</span>
        <span class="mf">2.2</span><span class="p">,</span> <span class="mf">2.4</span><span class="p">,</span> <span class="mf">2.6</span><span class="p">,</span> <span class="mf">2.8</span><span class="p">,</span> <span class="mf">3.</span><span class="p">,</span>
        <span class="mf">3.2</span><span class="p">,</span> <span class="mf">3.4</span><span class="p">,</span> <span class="mf">3.6</span><span class="p">,</span> <span class="mf">3.8</span>
        <span class="p">])</span>
    <span class="n">y</span> <span class="o">=</span> <span class="n">gv</span><span class="o">.</span><span class="n">gvar</span><span class="p">([</span>
        <span class="s1">&#39;0.38(20)&#39;</span><span class="p">,</span> <span class="s1">&#39;2.89(20)&#39;</span><span class="p">,</span> <span class="s1">&#39;0.85(20)&#39;</span><span class="p">,</span> <span class="s1">&#39;0.59(20)&#39;</span><span class="p">,</span> <span class="s1">&#39;2.88(20)&#39;</span><span class="p">,</span>
        <span class="s1">&#39;1.44(20)&#39;</span><span class="p">,</span> <span class="s1">&#39;0.73(20)&#39;</span><span class="p">,</span> <span class="s1">&#39;1.23(20)&#39;</span><span class="p">,</span> <span class="s1">&#39;1.68(20)&#39;</span><span class="p">,</span> <span class="s1">&#39;1.36(20)&#39;</span><span class="p">,</span>
        <span class="s1">&#39;1.51(20)&#39;</span><span class="p">,</span> <span class="s1">&#39;1.73(20)&#39;</span><span class="p">,</span> <span class="s1">&#39;2.16(20)&#39;</span><span class="p">,</span> <span class="s1">&#39;1.85(20)&#39;</span><span class="p">,</span> <span class="s1">&#39;2.00(20)&#39;</span><span class="p">,</span>
        <span class="s1">&#39;2.11(20)&#39;</span><span class="p">,</span> <span class="s1">&#39;2.75(20)&#39;</span><span class="p">,</span> <span class="s1">&#39;0.86(20)&#39;</span><span class="p">,</span> <span class="s1">&#39;2.73(20)&#39;</span>
        <span class="p">])</span>
    <span class="n">prior</span> <span class="o">=</span> <span class="n">make_prior</span><span class="p">()</span>
    <span class="n">fit</span> <span class="o">=</span> <span class="n">lsqfit</span><span class="o">.</span><span class="n">nonlinear_fit</span><span class="p">(</span><span class="n">data</span><span class="o">=</span><span class="p">(</span><span class="n">x</span><span class="p">,</span> <span class="n">y</span><span class="p">),</span> <span class="n">prior</span><span class="o">=</span><span class="n">prior</span><span class="p">,</span> <span class="n">fcn</span><span class="o">=</span><span class="n">fitfcn</span><span class="p">)</span>
    <span class="nb">print</span><span class="p">(</span><span class="n">fit</span><span class="p">)</span>

    <span class="c1">### 2) Bayesian integral with modified PDF</span>
    <span class="c1"># modified probability density function</span>
    <span class="n">mod_pdf</span> <span class="o">=</span> <span class="n">ModifiedPDF</span><span class="p">(</span><span class="n">data</span><span class="o">=</span><span class="p">(</span><span class="n">x</span><span class="p">,</span> <span class="n">y</span><span class="p">),</span> <span class="n">fcn</span><span class="o">=</span><span class="n">fitfcn</span><span class="p">,</span> <span class="n">prior</span><span class="o">=</span><span class="n">prior</span><span class="p">)</span>

    <span class="c1"># integrator for expectation values with modified PDF</span>
    <span class="n">expval</span> <span class="o">=</span> <span class="n">vegas</span><span class="o">.</span><span class="n">PDFIntegrator</span><span class="p">(</span><span class="n">fit</span><span class="o">.</span><span class="n">p</span><span class="p">,</span> <span class="n">pdf</span><span class="o">=</span><span class="n">mod_pdf</span><span class="p">)</span>

    <span class="c1"># adapt integrator to pdf</span>
    <span class="n">expval</span><span class="p">(</span><span class="n">neval</span><span class="o">=</span><span class="mi">1000</span><span class="p">,</span> <span class="n">nitn</span><span class="o">=</span><span class="mi">15</span><span class="p">)</span>

    <span class="c1"># evaluate expectation value of g(p)</span>
    <span class="k">def</span> <span class="nf">g</span><span class="p">(</span><span class="n">p</span><span class="p">):</span>
        <span class="n">w</span> <span class="o">=</span> <span class="n">p</span><span class="p">[</span><span class="s1">&#39;w&#39;</span><span class="p">]</span>
        <span class="n">c</span> <span class="o">=</span> <span class="n">p</span><span class="p">[</span><span class="s1">&#39;c&#39;</span><span class="p">]</span>
        <span class="k">return</span> <span class="nb">dict</span><span class="p">(</span><span class="n">w</span><span class="o">=</span><span class="p">[</span><span class="n">w</span><span class="p">,</span> <span class="n">w</span><span class="o">**</span><span class="mi">2</span><span class="p">],</span> <span class="n">mean</span><span class="o">=</span><span class="n">c</span><span class="p">,</span> <span class="n">outer</span><span class="o">=</span><span class="n">np</span><span class="o">.</span><span class="n">outer</span><span class="p">(</span><span class="n">c</span><span class="p">,</span><span class="n">c</span><span class="p">))</span>

    <span class="n">results</span> <span class="o">=</span> <span class="n">expval</span><span class="p">(</span><span class="n">g</span><span class="p">,</span> <span class="n">neval</span><span class="o">=</span><span class="mi">1000</span><span class="p">,</span> <span class="n">nitn</span><span class="o">=</span><span class="mi">15</span><span class="p">,</span> <span class="n">adapt</span><span class="o">=</span><span class="kc">False</span><span class="p">)</span>
    <span class="nb">print</span><span class="p">(</span><span class="n">results</span><span class="o">.</span><span class="n">summary</span><span class="p">())</span>

    <span class="c1"># parameters c[i]</span>
    <span class="n">mean</span> <span class="o">=</span> <span class="n">results</span><span class="p">[</span><span class="s1">&#39;mean&#39;</span><span class="p">]</span>
    <span class="n">cov</span> <span class="o">=</span> <span class="n">results</span><span class="p">[</span><span class="s1">&#39;outer&#39;</span><span class="p">]</span> <span class="o">-</span> <span class="n">np</span><span class="o">.</span><span class="n">outer</span><span class="p">(</span><span class="n">mean</span><span class="p">,</span> <span class="n">mean</span><span class="p">)</span>
    <span class="n">c</span> <span class="o">=</span> <span class="n">mean</span> <span class="o">+</span> <span class="n">gv</span><span class="o">.</span><span class="n">gvar</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">(</span><span class="n">mean</span><span class="o">.</span><span class="n">shape</span><span class="p">),</span> <span class="n">gv</span><span class="o">.</span><span class="n">mean</span><span class="p">(</span><span class="n">cov</span><span class="p">))</span>
    <span class="nb">print</span><span class="p">(</span><span class="s1">&#39;c =&#39;</span><span class="p">,</span> <span class="n">c</span><span class="p">)</span>
    <span class="nb">print</span><span class="p">(</span><span class="s1">&#39;corr(c) =&#39;</span><span class="p">,</span> <span class="nb">str</span><span class="p">(</span><span class="n">gv</span><span class="o">.</span><span class="n">evalcorr</span><span class="p">(</span><span class="n">c</span><span class="p">))</span><span class="o">.</span><span class="n">replace</span><span class="p">(</span><span class="s1">&#39;</span><span class="se">\n</span><span class="s1">&#39;</span><span class="p">,</span> <span class="s1">&#39;</span><span class="se">\n</span><span class="s1">&#39;</span> <span class="o">+</span> <span class="mi">10</span><span class="o">*</span><span class="s1">&#39; &#39;</span><span class="p">))</span>
    <span class="nb">print</span><span class="p">()</span>

    <span class="c1"># parameter w</span>
    <span class="n">wmean</span><span class="p">,</span> <span class="n">w2mean</span> <span class="o">=</span> <span class="n">results</span><span class="p">[</span><span class="s1">&#39;w&#39;</span><span class="p">]</span>
    <span class="n">wsdev</span> <span class="o">=</span> <span class="n">gv</span><span class="o">.</span><span class="n">mean</span><span class="p">(</span><span class="n">w2mean</span> <span class="o">-</span> <span class="n">wmean</span> <span class="o">**</span> <span class="mi">2</span><span class="p">)</span> <span class="o">**</span> <span class="mf">0.5</span>
    <span class="n">w</span> <span class="o">=</span> <span class="n">wmean</span> <span class="o">+</span> <span class="n">gv</span><span class="o">.</span><span class="n">gvar</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">shape</span><span class="p">(</span><span class="n">wmean</span><span class="p">)),</span> <span class="n">wsdev</span><span class="p">)</span>
    <span class="nb">print</span><span class="p">(</span><span class="s1">&#39;w =&#39;</span><span class="p">,</span> <span class="n">w</span><span class="p">,</span> <span class="s1">&#39;</span><span class="se">\n</span><span class="s1">&#39;</span><span class="p">)</span>

    <span class="c1"># Bayes Factor</span>
    <span class="nb">print</span><span class="p">(</span><span class="s1">&#39;logBF =&#39;</span><span class="p">,</span> <span class="n">np</span><span class="o">.</span><span class="n">log</span><span class="p">(</span><span class="n">results</span><span class="o">.</span><span class="n">pdfnorm</span><span class="p">))</span>

<span class="k">def</span> <span class="nf">fitfcn</span><span class="p">(</span><span class="n">x</span><span class="p">,</span> <span class="n">p</span><span class="p">):</span>
    <span class="n">c</span> <span class="o">=</span> <span class="n">p</span><span class="p">[</span><span class="s1">&#39;c&#39;</span><span class="p">]</span>
    <span class="k">return</span> <span class="n">c</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span> <span class="o">+</span> <span class="n">c</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span> <span class="o">*</span> <span class="n">x</span>

<span class="k">def</span> <span class="nf">make_prior</span><span class="p">():</span>
    <span class="n">prior</span> <span class="o">=</span> <span class="n">gv</span><span class="o">.</span><span class="n">BufferDict</span><span class="p">(</span><span class="n">c</span><span class="o">=</span><span class="n">gv</span><span class="o">.</span><span class="n">gvar</span><span class="p">([</span><span class="s1">&#39;0(5)&#39;</span><span class="p">,</span> <span class="s1">&#39;0(5)&#39;</span><span class="p">]))</span>
    <span class="n">prior</span><span class="p">[</span><span class="s1">&#39;gw(w)&#39;</span><span class="p">]</span> <span class="o">=</span> <span class="n">gv</span><span class="o">.</span><span class="n">BufferDict</span><span class="o">.</span><span class="n">uniform</span><span class="p">(</span><span class="s1">&#39;gw&#39;</span><span class="p">,</span> <span class="mf">0.</span><span class="p">,</span> <span class="mf">1.</span><span class="p">)</span>
    <span class="k">return</span> <span class="n">prior</span>

<span class="k">class</span> <span class="nc">ModifiedPDF</span><span class="p">:</span>
    <span class="sd">&quot;&quot;&quot; Modified PDF to account for measurement failure. &quot;&quot;&quot;</span>

    <span class="k">def</span> <span class="fm">__init__</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">data</span><span class="p">,</span> <span class="n">fcn</span><span class="p">,</span> <span class="n">prior</span><span class="p">):</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">x</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">y</span> <span class="o">=</span> <span class="n">data</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">fcn</span> <span class="o">=</span> <span class="n">fcn</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">prior</span> <span class="o">=</span> <span class="n">prior</span>

    <span class="k">def</span> <span class="fm">__call__</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">p</span><span class="p">):</span>
        <span class="n">w</span> <span class="o">=</span> <span class="n">p</span><span class="p">[</span><span class="s1">&#39;w&#39;</span><span class="p">]</span>
        <span class="n">y_fx</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">y</span> <span class="o">-</span> <span class="bp">self</span><span class="o">.</span><span class="n">fcn</span><span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">x</span><span class="p">,</span> <span class="n">p</span><span class="p">)</span>
        <span class="n">data_pdf1</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">gaussian_pdf</span><span class="p">(</span><span class="n">y_fx</span><span class="p">,</span> <span class="mf">1.</span><span class="p">)</span>
        <span class="n">data_pdf2</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">gaussian_pdf</span><span class="p">(</span><span class="n">y_fx</span><span class="p">,</span> <span class="mf">10.</span><span class="p">)</span>
        <span class="n">prior_pdf</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">gaussian_pdf</span><span class="p">(</span>
            <span class="n">p</span><span class="o">.</span><span class="n">buf</span><span class="p">[:</span><span class="nb">len</span><span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">prior</span><span class="o">.</span><span class="n">buf</span><span class="p">)]</span> <span class="o">-</span> <span class="bp">self</span><span class="o">.</span><span class="n">prior</span><span class="o">.</span><span class="n">buf</span>
            <span class="p">)</span>
        <span class="k">return</span> <span class="n">np</span><span class="o">.</span><span class="n">prod</span><span class="p">((</span><span class="mf">1.</span> <span class="o">-</span> <span class="n">w</span><span class="p">)</span> <span class="o">*</span> <span class="n">data_pdf1</span> <span class="o">+</span> <span class="n">w</span> <span class="o">*</span> <span class="n">data_pdf2</span><span class="p">)</span> <span class="o">*</span> <span class="n">np</span><span class="o">.</span><span class="n">prod</span><span class="p">(</span><span class="n">prior_pdf</span><span class="p">)</span>

    <span class="nd">@staticmethod</span>
    <span class="k">def</span> <span class="nf">gaussian_pdf</span><span class="p">(</span><span class="n">x</span><span class="p">,</span> <span class="n">f</span><span class="o">=</span><span class="mf">1.</span><span class="p">):</span>
        <span class="n">xmean</span> <span class="o">=</span> <span class="n">gv</span><span class="o">.</span><span class="n">mean</span><span class="p">(</span><span class="n">x</span><span class="p">)</span>
        <span class="n">xvar</span> <span class="o">=</span> <span class="n">gv</span><span class="o">.</span><span class="n">var</span><span class="p">(</span><span class="n">x</span><span class="p">)</span> <span class="o">*</span> <span class="n">f</span> <span class="o">**</span> <span class="mi">2</span>
        <span class="k">return</span> <span class="n">gv</span><span class="o">.</span><span class="n">exp</span><span class="p">(</span><span class="o">-</span><span class="n">xmean</span> <span class="o">**</span> <span class="mi">2</span> <span class="o">/</span> <span class="mf">2.</span> <span class="o">/</span><span class="n">xvar</span><span class="p">)</span> <span class="o">/</span> <span class="n">gv</span><span class="o">.</span><span class="n">sqrt</span><span class="p">(</span><span class="mi">2</span> <span class="o">*</span> <span class="n">np</span><span class="o">.</span><span class="n">pi</span> <span class="o">*</span> <span class="n">xvar</span><span class="p">)</span>

<span class="k">if</span> <span class="vm">__name__</span> <span class="o">==</span> <span class="s1">&#39;__main__&#39;</span><span class="p">:</span>
    <span class="n">main</span><span class="p">()</span>
</pre></div>
</div>
<p>Here class <code class="docutils literal notranslate"><span class="pre">ModifiedPDF</span></code> implements the modified PDF.  As usual the PDF for
the parameters (in <code class="docutils literal notranslate"><span class="pre">__call__</span></code>) is the product of a PDF for the data times a
PDF for the priors. The data PDF is more complicated than usual, however, as
it consists of two Gaussian distributions: one, <code class="docutils literal notranslate"><span class="pre">data_pdf1</span></code>, with the
nominal data errors, and the other, <code class="docutils literal notranslate"><span class="pre">data_pdf2</span></code>, with errors that are ten
times larger. Parameter <code class="docutils literal notranslate"><span class="pre">w</span></code> determines the relative weight of each data PDF.</p>
<p>The Bayesian integrals are estimated using <code class="xref py py-class docutils literal notranslate"><span class="pre">vegas.PDFIntegrator</span></code>
<code class="docutils literal notranslate"><span class="pre">expval</span></code>, which is created from the least-squares fit output (<code class="docutils literal notranslate"><span class="pre">fit</span></code>).
It is used to evaluate expectation values of arbitrary functions of the
fit variables. Normally it would use the standard PDF from the least-squares
fit, but we replace that PDF here with an instance (<code class="docutils literal notranslate"><span class="pre">mod_pdf</span></code>) of class
<code class="docutils literal notranslate"><span class="pre">ModifiedPDF</span></code>.</p>
<p>We have modified <code class="docutils literal notranslate"><span class="pre">make_prior()</span></code> to introduce <code class="docutils literal notranslate"><span class="pre">w</span></code> as a new fit
parameter. The prior for this parameter is uniformly distributed
across the interval from 0 to 1. Parameter <code class="docutils literal notranslate"><span class="pre">w</span></code> plays
no role in the initial fit. (The uniform distribution is implemented
by introducing a function <code class="docutils literal notranslate"><span class="pre">gw(w)</span></code> that maps it onto a Gaussian
distribution 0 ± 1. The integration parameter in the
Bayesian integrals is <code class="docutils literal notranslate"><span class="pre">gw(w)</span></code> but the <code class="docutils literal notranslate"><span class="pre">BufferDict</span></code> dictionary
makes the corresponding value of <code class="docutils literal notranslate"><span class="pre">w</span></code> available automatically.)</p>
<p>We first call <code class="docutils literal notranslate"><span class="pre">expval</span></code> with no function, to allow the integrator to adapt
to the modified PDF. We then use the integrator, now with adaptation
turned off (<code class="docutils literal notranslate"><span class="pre">adapt=False</span></code>), to evaluate the expectation value of
function <code class="docutils literal notranslate"><span class="pre">g(p)</span></code>. The output dictionary <code class="docutils literal notranslate"><span class="pre">results</span></code>
contains expectation values of the corresponding entries in the dictionary
returned <code class="docutils literal notranslate"><span class="pre">g(p)</span></code>. These data allow us to calculate means, standard deviations
and correlation matrices for the fit parameters.</p>
<p>The results from this code are as follows:</p>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span><span class="n">Least</span> <span class="n">Square</span> <span class="n">Fit</span><span class="p">:</span>
  <span class="n">chi2</span><span class="o">/</span><span class="n">dof</span> <span class="p">[</span><span class="n">dof</span><span class="p">]</span> <span class="o">=</span> <span class="mi">13</span> <span class="p">[</span><span class="mi">19</span><span class="p">]</span>    <span class="n">Q</span> <span class="o">=</span> <span class="mf">1.2e-40</span>    <span class="n">logGBF</span> <span class="o">=</span> <span class="o">-</span><span class="mf">117.45</span>

<span class="n">Parameters</span><span class="p">:</span>
            <span class="n">c</span> <span class="mi">0</span>   <span class="mf">1.149</span> <span class="p">(</span><span class="mi">95</span><span class="p">)</span>     <span class="p">[</span>  <span class="mf">0.0</span> <span class="p">(</span><span class="mf">5.0</span><span class="p">)</span> <span class="p">]</span>  
              <span class="mi">1</span>   <span class="mf">0.261</span> <span class="p">(</span><span class="mi">42</span><span class="p">)</span>     <span class="p">[</span>  <span class="mf">0.0</span> <span class="p">(</span><span class="mf">5.0</span><span class="p">)</span> <span class="p">]</span>  
          <span class="n">gw</span><span class="p">(</span><span class="n">w</span><span class="p">)</span>   <span class="mf">8e-15</span> <span class="o">+-</span> <span class="mi">1</span>     <span class="p">[</span>  <span class="mf">0.0</span> <span class="p">(</span><span class="mf">1.0</span><span class="p">)</span> <span class="p">]</span>  
<span class="o">-----------------------------------------------</span>
              <span class="n">w</span>    <span class="mf">0.50</span> <span class="p">(</span><span class="mi">40</span><span class="p">)</span>     <span class="p">[</span>  <span class="mf">0.50</span> <span class="p">(</span><span class="mi">40</span><span class="p">)</span> <span class="p">]</span>  

<span class="n">Settings</span><span class="p">:</span>
  <span class="n">svdcut</span><span class="o">/</span><span class="n">n</span> <span class="o">=</span> <span class="mf">1e-12</span><span class="o">/</span><span class="mi">0</span>    <span class="n">tol</span> <span class="o">=</span> <span class="p">(</span><span class="mf">1e-08</span><span class="o">*</span><span class="p">,</span><span class="mf">1e-10</span><span class="p">,</span><span class="mf">1e-10</span><span class="p">)</span>    <span class="p">(</span><span class="n">itns</span><span class="o">/</span><span class="n">time</span> <span class="o">=</span> <span class="mi">5</span><span class="o">/</span><span class="mf">0.0</span><span class="p">)</span>

<span class="n">itn</span>   <span class="n">integral</span>        <span class="n">average</span>         <span class="n">chi2</span><span class="o">/</span><span class="n">dof</span>        <span class="n">Q</span>
<span class="o">-------------------------------------------------------</span>
  <span class="mi">1</span>   <span class="mf">6.821</span><span class="p">(</span><span class="mi">73</span><span class="p">)</span><span class="n">e</span><span class="o">-</span><span class="mi">11</span>   <span class="mf">6.821</span><span class="p">(</span><span class="mi">73</span><span class="p">)</span><span class="n">e</span><span class="o">-</span><span class="mi">11</span>       <span class="mf">0.00</span>     <span class="mf">1.00</span>
  <span class="mi">2</span>   <span class="mf">6.849</span><span class="p">(</span><span class="mi">65</span><span class="p">)</span><span class="n">e</span><span class="o">-</span><span class="mi">11</span>   <span class="mf">6.835</span><span class="p">(</span><span class="mi">49</span><span class="p">)</span><span class="n">e</span><span class="o">-</span><span class="mi">11</span>       <span class="mf">0.62</span>     <span class="mf">0.78</span>
  <span class="mi">3</span>   <span class="mf">6.855</span><span class="p">(</span><span class="mi">75</span><span class="p">)</span><span class="n">e</span><span class="o">-</span><span class="mi">11</span>   <span class="mf">6.842</span><span class="p">(</span><span class="mi">41</span><span class="p">)</span><span class="n">e</span><span class="o">-</span><span class="mi">11</span>       <span class="mf">0.62</span>     <span class="mf">0.89</span>
  <span class="mi">4</span>   <span class="mf">6.853</span><span class="p">(</span><span class="mi">64</span><span class="p">)</span><span class="n">e</span><span class="o">-</span><span class="mi">11</span>   <span class="mf">6.845</span><span class="p">(</span><span class="mi">35</span><span class="p">)</span><span class="n">e</span><span class="o">-</span><span class="mi">11</span>       <span class="mf">0.51</span>     <span class="mf">0.98</span>
  <span class="mi">5</span>   <span class="mf">6.731</span><span class="p">(</span><span class="mi">70</span><span class="p">)</span><span class="n">e</span><span class="o">-</span><span class="mi">11</span>   <span class="mf">6.822</span><span class="p">(</span><span class="mi">31</span><span class="p">)</span><span class="n">e</span><span class="o">-</span><span class="mi">11</span>       <span class="mf">0.56</span>     <span class="mf">0.99</span>
  <span class="mi">6</span>   <span class="mf">6.827</span><span class="p">(</span><span class="mi">83</span><span class="p">)</span><span class="n">e</span><span class="o">-</span><span class="mi">11</span>   <span class="mf">6.823</span><span class="p">(</span><span class="mi">29</span><span class="p">)</span><span class="n">e</span><span class="o">-</span><span class="mi">11</span>       <span class="mf">0.77</span>     <span class="mf">0.87</span>
  <span class="mi">7</span>   <span class="mf">6.734</span><span class="p">(</span><span class="mi">78</span><span class="p">)</span><span class="n">e</span><span class="o">-</span><span class="mi">11</span>   <span class="mf">6.810</span><span class="p">(</span><span class="mi">28</span><span class="p">)</span><span class="n">e</span><span class="o">-</span><span class="mi">11</span>       <span class="mf">0.77</span>     <span class="mf">0.89</span>
  <span class="mi">8</span>   <span class="mf">6.752</span><span class="p">(</span><span class="mi">79</span><span class="p">)</span><span class="n">e</span><span class="o">-</span><span class="mi">11</span>   <span class="mf">6.803</span><span class="p">(</span><span class="mi">26</span><span class="p">)</span><span class="n">e</span><span class="o">-</span><span class="mi">11</span>       <span class="mf">0.77</span>     <span class="mf">0.91</span>
  <span class="mi">9</span>   <span class="mf">6.816</span><span class="p">(</span><span class="mi">73</span><span class="p">)</span><span class="n">e</span><span class="o">-</span><span class="mi">11</span>   <span class="mf">6.804</span><span class="p">(</span><span class="mi">25</span><span class="p">)</span><span class="n">e</span><span class="o">-</span><span class="mi">11</span>       <span class="mf">0.71</span>     <span class="mf">0.97</span>
 <span class="mi">10</span>   <span class="mf">6.740</span><span class="p">(</span><span class="mi">76</span><span class="p">)</span><span class="n">e</span><span class="o">-</span><span class="mi">11</span>   <span class="mf">6.798</span><span class="p">(</span><span class="mi">23</span><span class="p">)</span><span class="n">e</span><span class="o">-</span><span class="mi">11</span>       <span class="mf">0.77</span>     <span class="mf">0.94</span>
 <span class="mi">11</span>   <span class="mf">6.843</span><span class="p">(</span><span class="mi">66</span><span class="p">)</span><span class="n">e</span><span class="o">-</span><span class="mi">11</span>   <span class="mf">6.802</span><span class="p">(</span><span class="mi">22</span><span class="p">)</span><span class="n">e</span><span class="o">-</span><span class="mi">11</span>       <span class="mf">0.81</span>     <span class="mf">0.91</span>
 <span class="mi">12</span>   <span class="mf">6.729</span><span class="p">(</span><span class="mi">67</span><span class="p">)</span><span class="n">e</span><span class="o">-</span><span class="mi">11</span>   <span class="mf">6.796</span><span class="p">(</span><span class="mi">21</span><span class="p">)</span><span class="n">e</span><span class="o">-</span><span class="mi">11</span>       <span class="mf">0.81</span>     <span class="mf">0.91</span>
 <span class="mi">13</span>   <span class="mf">6.746</span><span class="p">(</span><span class="mi">71</span><span class="p">)</span><span class="n">e</span><span class="o">-</span><span class="mi">11</span>   <span class="mf">6.792</span><span class="p">(</span><span class="mi">20</span><span class="p">)</span><span class="n">e</span><span class="o">-</span><span class="mi">11</span>       <span class="mf">0.78</span>     <span class="mf">0.95</span>
 <span class="mi">14</span>   <span class="mf">6.712</span><span class="p">(</span><span class="mi">78</span><span class="p">)</span><span class="n">e</span><span class="o">-</span><span class="mi">11</span>   <span class="mf">6.786</span><span class="p">(</span><span class="mi">19</span><span class="p">)</span><span class="n">e</span><span class="o">-</span><span class="mi">11</span>       <span class="mf">0.78</span>     <span class="mf">0.96</span>
 <span class="mi">15</span>   <span class="mf">6.766</span><span class="p">(</span><span class="mi">72</span><span class="p">)</span><span class="n">e</span><span class="o">-</span><span class="mi">11</span>   <span class="mf">6.785</span><span class="p">(</span><span class="mi">19</span><span class="p">)</span><span class="n">e</span><span class="o">-</span><span class="mi">11</span>       <span class="mf">0.79</span>     <span class="mf">0.96</span>

<span class="n">c</span> <span class="o">=</span> <span class="p">[</span><span class="mf">0.28</span><span class="p">(</span><span class="mi">14</span><span class="p">)</span> <span class="mf">0.623</span><span class="p">(</span><span class="mi">58</span><span class="p">)]</span>
<span class="n">corr</span><span class="p">(</span><span class="n">c</span><span class="p">)</span> <span class="o">=</span> <span class="p">[[</span> <span class="mf">1.</span>        <span class="o">-</span><span class="mf">0.8981073</span><span class="p">]</span>
           <span class="p">[</span><span class="o">-</span><span class="mf">0.8981073</span>  <span class="mf">1.</span>       <span class="p">]]</span>

<span class="n">w</span> <span class="o">=</span> <span class="mf">0.26</span><span class="p">(</span><span class="mi">11</span><span class="p">)</span> 

<span class="n">logBF</span> <span class="o">=</span> <span class="o">-</span><span class="mf">23.4137</span><span class="p">(</span><span class="mi">28</span><span class="p">)</span>
</pre></div>
</div>
<p>The table after the fit shows results for the normalization of the
modified PDF from each of the <code class="docutils literal notranslate"><span class="pre">nitn=15</span></code> iterations of the <code class="xref py py-mod docutils literal notranslate"><span class="pre">vegas</span></code>
algorithm used to estimate the integrals. The logarithm of the normalization
(<code class="docutils literal notranslate"><span class="pre">logBF</span></code>) is -23.4, which is much larger than the value -117.5 of <code class="docutils literal notranslate"><span class="pre">logGBF</span></code>
from the least-squares fit. This means that the data much prefer the
modified prior (by a factor of <code class="docutils literal notranslate"><span class="pre">exp(-23.4</span> <span class="pre">+</span> <span class="pre">117.4)</span></code> or about 10<sup>41</sup>.).</p>
<p>The new fit parameters are much more reasonable. In particular the
intercept is 0.28(14) rather than the 1.15(10) from the least-squares fit.
This is much better suited to the data (see the dashed line in red, with
the red band showing the 1-sigma region about the best fit):</p>
<a class="reference internal image-reference" href="_images/case-outliers2.png"><img alt="_images/case-outliers2.png" src="_images/case-outliers2.png" style="width: 60%;" /></a>
<p>Note, from the correlation matrix, that the intercept and slope are
anti-correlated, as one might guess for this fit.
The analysis also gives us an estimate for the failure rate <code class="docutils literal notranslate"><span class="pre">w=0.26(11)</span></code>
of our detectors — they fail about a quarter of the time.</p>
</section>
<section id="a-variation">
<h2>A Variation<a class="headerlink" href="#a-variation" title="Permalink to this heading">¶</a></h2>
<p>A slightly different model for the failure that
leads to outliers assigns a different <code class="docutils literal notranslate"><span class="pre">w</span></code> to each data point.
It is easily implemented here by changing the prior
so that <code class="docutils literal notranslate"><span class="pre">w</span></code> is an array:</p>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span><span class="k">def</span> <span class="nf">make_prior</span><span class="p">():</span>
    <span class="n">prior</span> <span class="o">=</span> <span class="n">gv</span><span class="o">.</span><span class="n">BufferDict</span><span class="p">(</span><span class="n">c</span><span class="o">=</span><span class="n">gv</span><span class="o">.</span><span class="n">gvar</span><span class="p">([</span><span class="s1">&#39;0(5)&#39;</span><span class="p">,</span> <span class="s1">&#39;0(5)&#39;</span><span class="p">]))</span>
    <span class="n">prior</span><span class="p">[</span><span class="s1">&#39;gw(w)&#39;</span><span class="p">]</span> <span class="o">=</span> <span class="n">gv</span><span class="o">.</span><span class="n">BufferDict</span><span class="o">.</span><span class="n">uniform</span><span class="p">(</span><span class="s1">&#39;gw&#39;</span><span class="p">,</span> <span class="mf">0.</span><span class="p">,</span> <span class="mf">1.</span><span class="p">,</span> <span class="n">shape</span><span class="o">=</span><span class="mi">19</span><span class="p">)</span>
    <span class="k">return</span> <span class="n">prior</span>
</pre></div>
</div>
<p>The Bayesian integral then has 21 parameters, rather than the 3 parameters
before. The code still takes only 4 secs to run (on a 2020 laptop).</p>
<p>The final results are quite similar to the other model:</p>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span><span class="n">c</span> <span class="o">=</span> <span class="p">[</span><span class="mf">0.29</span><span class="p">(</span><span class="mi">16</span><span class="p">)</span> <span class="mf">0.616</span><span class="p">(</span><span class="mi">67</span><span class="p">)]</span>
<span class="n">corr</span><span class="p">(</span><span class="n">c</span><span class="p">)</span> <span class="o">=</span> <span class="p">[[</span> <span class="mf">1.</span>         <span class="o">-</span><span class="mf">0.89773341</span><span class="p">]</span>
           <span class="p">[</span><span class="o">-</span><span class="mf">0.89773341</span>  <span class="mf">1.</span>        <span class="p">]]</span> 

<span class="n">w</span> <span class="o">=</span> <span class="p">[</span><span class="mf">0.37</span><span class="p">(</span><span class="mi">26</span><span class="p">)</span> <span class="mf">0.67</span><span class="p">(</span><span class="mi">23</span><span class="p">)</span> <span class="mf">0.39</span><span class="p">(</span><span class="mi">27</span><span class="p">)</span> <span class="mf">0.38</span><span class="p">(</span><span class="mi">26</span><span class="p">)</span> <span class="mf">0.68</span><span class="p">(</span><span class="mi">23</span><span class="p">)</span> <span class="mf">0.49</span><span class="p">(</span><span class="mi">29</span><span class="p">)</span> <span class="mf">0.50</span><span class="p">(</span><span class="mi">29</span><span class="p">)</span> <span class="mf">0.37</span><span class="p">(</span><span class="mi">25</span><span class="p">)</span>
 <span class="mf">0.43</span><span class="p">(</span><span class="mi">28</span><span class="p">)</span> <span class="mf">0.37</span><span class="p">(</span><span class="mi">26</span><span class="p">)</span> <span class="mf">0.37</span><span class="p">(</span><span class="mi">27</span><span class="p">)</span> <span class="mf">0.37</span><span class="p">(</span><span class="mi">26</span><span class="p">)</span> <span class="mf">0.42</span><span class="p">(</span><span class="mi">28</span><span class="p">)</span> <span class="mf">0.38</span><span class="p">(</span><span class="mi">26</span><span class="p">)</span> <span class="mf">0.38</span><span class="p">(</span><span class="mi">26</span><span class="p">)</span> <span class="mf">0.39</span><span class="p">(</span><span class="mi">27</span><span class="p">)</span>
 <span class="mf">0.47</span><span class="p">(</span><span class="mi">29</span><span class="p">)</span> <span class="mf">0.67</span><span class="p">(</span><span class="mi">24</span><span class="p">)</span> <span class="mf">0.38</span><span class="p">(</span><span class="mi">26</span><span class="p">)]</span> 

<span class="n">logBF</span> <span class="o">=</span> <span class="o">-</span><span class="mf">24.125</span><span class="p">(</span><span class="mi">23</span><span class="p">)</span>
</pre></div>
</div>
<p>Note that the logarithm of the Bayes Factor <code class="docutils literal notranslate"><span class="pre">logBF</span></code> is slighly lower for
this model than before. It is also less accurately determined (20x), because
21-parameter integrals are considerably more difficult than 3-parameter
integrals. More precision can be obtained by increasing <code class="docutils literal notranslate"><span class="pre">neval</span></code>, but
the current precision is more than adequate.</p>
<p>Only three of the <code class="docutils literal notranslate"><span class="pre">w[i]</span></code> values listed in the output are more than two
standard deviations away from zero. Not surprisingly, these correspond to
the unambiguous outliers.</p>
<p>The outliers in this case are pretty obvious; one is tempted to simply drop
them. It is clearly better, however, to understand why they have occurred and
to quantify the effect if possible, as above. Dropping outliers would be much
more difficult if they were, say, three times closer to the rest of the data.
The least-squares fit would still be poor (<code class="docutils literal notranslate"><span class="pre">chi**2</span></code> per degree of freedom of
3) and its intercept a bit too high (0.6(1)). Using the modified PDF, on the
other hand, would give results very similar to what we obtained above: for
example, the intercept would be 0.35(17).</p>
</section>
</section>


            <div class="clearer"></div>
          </div>
        </div>
      </div>
      <div class="sphinxsidebar" role="navigation" aria-label="main navigation">
        <div class="sphinxsidebarwrapper">
  <div>
    <h3><a href="index.html">Table of Contents</a></h3>
    <ul>
<li><a class="reference internal" href="#">Case Study: Outliers and Bayesian Integrals</a><ul>
<li><a class="reference internal" href="#the-problem">The Problem</a></li>
<li><a class="reference internal" href="#a-solution">A Solution</a></li>
<li><a class="reference internal" href="#a-variation">A Variation</a></li>
</ul>
</li>
</ul>

  </div>
  <div>
    <h4>Previous topic</h4>
    <p class="topless"><a href="case-spline.html"
                          title="previous chapter">Case Study: Fitting a Spline</a></p>
  </div>
  <div>
    <h4>Next topic</h4>
    <p class="topless"><a href="lsqfit.html"
                          title="next chapter"><code class="xref py py-mod docutils literal notranslate"><span class="pre">lsqfit</span></code> - Nonlinear Least Squares Fitting</a></p>
  </div>
  <div role="note" aria-label="source link">
    <h3>This Page</h3>
    <ul class="this-page-menu">
      <li><a href="_sources/case-outliers.rst.txt"
            rel="nofollow">Show Source</a></li>
    </ul>
   </div>
<div id="searchbox" style="display: none" role="search">
  <h3 id="searchlabel">Quick search</h3>
    <div class="searchformwrapper">
    <form class="search" action="search.html" method="get">
      <input type="text" name="q" aria-labelledby="searchlabel" autocomplete="off" autocorrect="off" autocapitalize="off" spellcheck="false"/>
      <input type="submit" value="Go" />
    </form>
    </div>
</div>
<script>document.getElementById('searchbox').style.display = "block"</script>
        </div>
      </div>
      <div class="clearer"></div>
    </div>
    <div class="related" role="navigation" aria-label="related navigation">
      <h3>Navigation</h3>
      <ul>
        <li class="right" style="margin-right: 10px">
          <a href="genindex.html" title="General Index"
             >index</a></li>
        <li class="right" >
          <a href="py-modindex.html" title="Python Module Index"
             >modules</a> |</li>
        <li class="right" >
          <a href="lsqfit.html" title="lsqfit - Nonlinear Least Squares Fitting"
             >next</a> |</li>
        <li class="right" >
          <a href="case-spline.html" title="Case Study: Fitting a Spline"
             >previous</a> |</li>
        <li class="nav-item nav-item-0"><a href="index.html">lsqfit 13.0.1 documentation</a> &#187;</li>
        <li class="nav-item nav-item-this"><a href="">Case Study: Outliers and Bayesian Integrals</a></li> 
      </ul>
    </div>
    <div class="footer" role="contentinfo">
        &#169; Copyright 2009-2019, G. P. Lepage.
      Created using <a href="https://www.sphinx-doc.org/">Sphinx</a> 5.0.2.
    </div>
  </body>
</html>