tag:blogger.com,1999:blog-13694323968982046132024-03-06T00:28:13.176+01:00Lighthouse in the SkyAstronomy, science, technology, whimsy. Not necessarily in that order.Unknownnoreply@blogger.comBlogger138125tag:blogger.com,1999:blog-1369432396898204613.post-58920703956441443682017-04-18T20:00:00.000+02:002017-04-20T20:02:43.656+02:00Curve Fitting Part 6: Linear Least Squares and its Errors <div tabindex="-1" id="notebook" class="border-box-sizing">
<div class="container" id="notebook-container">
<div class="cell border-box-sizing text_cell rendered">
<div class="prompt input_prompt">
</div>
<div class="inner_cell">
<div class="text_cell_render border-box-sizing rendered_html">
<h1 id="linear-least-squares-and-its-uncertainties">Linear least squares and its uncertainties</h1>
</div>
</div>
</div>
<div class="cell border-box-sizing code_cell rendered">
<div class="input">
<div class="prompt input_prompt">In [1]:</div>
<div class="inner_cell">
<div class="input_area">
<div class="highlight"><pre><span class="kn">from</span> <span class="nn">collections</span> <span class="kn">import</span> <span class="n">namedtuple</span>
<span class="kn">import</span> <span class="nn">numpy</span> <span class="kn">as</span> <span class="nn">np</span>
<span class="kn">import</span> <span class="nn">numpy.linalg</span>
<span class="kn">import</span> <span class="nn">matplotlib</span>
<span class="n">matplotlib</span><span class="o">.</span><span class="n">rcParams</span><span class="p">[</span><span class="s1">'savefig.dpi'</span><span class="p">]</span> <span class="o">=</span> <span class="mi">144</span>
<span class="n">matplotlib</span><span class="o">.</span><span class="n">rcParams</span><span class="p">[</span><span class="s2">"image.composite_image"</span><span class="p">]</span><span class="o">=</span><span class="bp">False</span>
<span class="o">%</span><span class="k">matplotlib</span> <span class="n">inline</span>
<span class="kn">import</span> <span class="nn">matplotlib.pyplot</span> <span class="kn">as</span> <span class="nn">plt</span>
</pre></div>
</div>
</div>
</div>
</div>
<div class="cell border-box-sizing text_cell rendered">
<div class="prompt input_prompt">
</div>
<div class="inner_cell">
<div class="text_cell_render border-box-sizing rendered_html">
<h2 id="setting-up-a-simple-test-problem">Setting up a simple test problem</h2>
<p>The goal is to demonstrate techniques of error estimation for linear least squares, so let's set up a problem about as simple as possible while still presenting the difficulties we might hope these techniques solve.</p>
<p>The problem will be a simple straight-line fit to data. I'm not even going to bother fooling with the best-fit values; the slope and intercept will be zero, and the data points are symmetrical about zero, so the slope/intercept parameterization should yield uncorrelated errors on the fitted parameters.</p>
<p>The goal of the fitting is to find <span class="math inline">\(x\)</span> that minimizes <span class="math display">\[
\left\|Ax-b\right\|^2_u,
\]</span> that is, the sum of squares of the residuals <span class="math inline">\(r\)</span> divided by (our estimate of) the uncertainties <span class="math inline">\(u\)</span>. We will be interested in what happens when our estimate of the uncertainties differs from the true uncertainties, so here I make the distinction.</p>
</div>
</div>
</div>
<div class="cell border-box-sizing code_cell rendered">
<div class="input">
<div class="prompt input_prompt">In [2]:</div>
<div class="inner_cell">
<div class="input_area">
<div class="highlight"><pre><span class="n">xs</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">linspace</span><span class="p">(</span><span class="o">-</span><span class="mi">1</span><span class="p">,</span><span class="mi">1</span><span class="p">,</span><span class="mi">100</span><span class="p">)</span>
<span class="n">true_uncerts</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">abs</span><span class="p">(</span><span class="n">xs</span><span class="p">)</span>
<span class="n">claimed_uncerts</span> <span class="o">=</span> <span class="mf">0.35</span><span class="o">*</span><span class="n">np</span><span class="o">.</span><span class="n">ones_like</span><span class="p">(</span><span class="n">xs</span><span class="p">)</span>
<span class="k">def</span> <span class="nf">make_fake_data</span><span class="p">(</span><span class="n">xs</span><span class="p">,</span> <span class="n">true_uncerts</span><span class="p">):</span>
<span class="k">return</span> <span class="n">true_uncerts</span><span class="o">*</span><span class="n">np</span><span class="o">.</span><span class="n">random</span><span class="o">.</span><span class="n">randn</span><span class="p">(</span><span class="nb">len</span><span class="p">(</span><span class="n">xs</span><span class="p">))</span>
<span class="n">ys</span> <span class="o">=</span> <span class="n">make_fake_data</span><span class="p">(</span><span class="n">xs</span><span class="p">,</span> <span class="n">true_uncerts</span><span class="p">)</span>
<span class="n">plt</span><span class="o">.</span><span class="n">errorbar</span><span class="p">(</span><span class="n">xs</span><span class="p">,</span><span class="n">ys</span><span class="p">,</span><span class="n">true_uncerts</span><span class="p">,</span> <span class="n">linestyle</span><span class="o">=</span><span class="s1">'none'</span><span class="p">)</span>
</pre></div>
</div>
</div>
</div>
<div class="output_wrapper">
<div class="output">
<div class="output_area"><div class="prompt output_prompt">Out[2]:</div>
<div class="output_text output_subarea output_pyout">
<pre>
<Container object of 3 artists>
</pre>
</div>
</div>
<div class="output_area"><div class="prompt"></div>
<div class="output_png output_subarea ">
<img src="data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAXcAAAEACAYAAABI5zaHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz
AAALEgAACxIB0t1+/AAAFMhJREFUeJzt3W2sHNddx/Hf340S8SASStSkitsYEVpEKTRGKuVBdE0V
4aBWNqhI6RuaFqkgJEoFlUhJpXuNEIJISAH6AhChShBRqPIijUtb4hBvKh5CoyauHVqnLmraOHFS
QbGqUFRFyZ8Xu9d37/XO7NkzZ2bOnPl+pCuv78w9c3Z35rf/PfNk7i4AQFn29N0BAEB6hDsAFIhw
B4ACEe4AUCDCHQAKRLgDQIEuadqAmV0m6TOSLp23d6+7H2naLgAgnqU4zt3MvtPdv2Vmr5D0L5Le
7+6fbdwwACBKkmEZd//W/OFlmlXvnBkFAD1KEu5mtsfMHpf0nKRj7v5oinYBAHFSVe4vu/v1kvZK
+gkz++EU7QIA4jTeobrI3b9pZsclHZT0hcVpZsZQDQBEcHdb928aV+5mdqWZXT5//B2SbpB0etm8
7s5Pop+NjY3e+1DKD68lr2fOP7FSVO6vlnSnme3R7MPi7939kwnaBQBEahzu7n5K0v4EfQEAJMIZ
qgM1mUz67kIxeC3T4vXMQ5KTmIIWZOZdLQsASmFm8j52qAIA8kO4A0CBCHcAKBDhDgAFItwBoECE
OwAUiHAHgAIR7gBQIMIdAApEuANAgQh3ACgQ4Q4ABSLcAaBAhDsAFIhwB4ACEe4AUCDCHQAKRLgD
QIEa3yAbALoync5+th5v3a51Mtl+jJlO76G6sTFbFm8EgKbMpDHcljn2HqrcIBvAIBHu9RhzB4AC
Ee4AUCDCHQAKRLgDQIEIdwAoUONwN7O9ZvaQmf2HmZ0ys/en6BgAIF7jQyHN7GpJV7v7CTP7bkmf
k3TI3U/vmo9DIQEkw6GQ9Rqfoeruz0l6bv74BTP7oqRrJJ2u/cOB4cw4AEOS9CQmM9snaSrpR9z9
hV3Tiqncx1IxADkby3bY+0lM8yGZeyX91u5gBwB0K8mFw8zsEs2C/W/d/eNV821ubl54PJlMNGE8
AwB2mE6nmm6NATeQZFjGzO6S9F/u/ts18zAsAyCZsWyHvV04zMx+WtJnJJ2S5POf33P3T++aj3AH
kMxYtkOuCtmhsaxUQM7Gsh32vkMVAJAPwh0ACkS4A0CBRnEPVc4uBTA2o9uhmmInzFh25AA5G8t2
yA5VAMAFhDsAFKjTMfetqw8w7g0A7eplzL3PsTLG3IEyjGU7ZMwdAHAB4Q4ABSLcAaBAoziJCUCZ
OEGxGjtUG7bBygX0Y/e2nMsO1tSZMKhL/pYU7qnbBhAm13BflCZv4sI9i2EZql8ASCu7yr3tT18q
d6AMVO71sqjcAZSDb+J5oHJP2EaOlQPQpza3CSr3elTuiVG1AMhB9pV7+sOKuqvcc6wkgC5RufdX
uWcf7qHTwvtBuANdySXc+/pGTbjvmI9wB0qRS7h31ac2lsWYOwAEGsO9JajcE7ZB5Q7slHvlPoyj
86jcO7Xskx/A8JR6hBuVe8M26qoAKneM3dAq99T9pXIHgLlSK+muUbk3bIPKHajWdBtYJw+o3HdK
cicmM7vDzJ43s5Mp2gMANJOkcjezn5H0gqS73P1HK+ahcgdGpqTKPWa4KMUZ972fxGRm10o6Srin
XRaQo9CQKincQ/sVOl94G4R7YD8JdyClNrfZdSrfAwcI90WdHi2zuXVwuKTpdKIJu74BrGEMR8xM
p1NNtz7BGiimcu/qa+LuNkqu3DkkDSG6rNxTzDuWyj1luO/TLNzfWDG9s2GZLod2Sg73RSU9F6RF
uMfNN4hhGTO7W9JE0veZ2dckbbj7R1O0DWDc+AYZp8iTmKjc0yvpuSCtLiv31PnQpHLf2Jg9jh0G
XpxWv6OY67knbSO0fcIdYzfWcK9qI7RPVdMu7lOPZ6gCAPJCuANAgQh3ACgQl/wF0AuOgmkXO1Qj
sEMV2NbmNssO1YFcfiB3dZXE1u8k6a1v3b7NHgDkiMo9oo2q+ajcMUZU7nlW7uxQBYACMSzTsq3h
mxJ2GJX0XIAcLNumUmFYJqKNqvnqvuLVtRd61MDu+fbtk556avvfNgM39LlgfBiWaZ5n9X3i8gOd
tVE1X2y4p+hf24Eb81w41G0cCPc8w51hGbRmMcTNtoMe8XL9wKwashuyVUfP5Y7KPaKNqvmo3Ov/
/vjxPINpqGLf89QfEKu2gSZXT8ylco/dthmW2TEf4V5quHfZ3yHr665iKdpIHYKEO+HeaRtV8/Ud
7iEVUqzdz2XdZRHucVKsy7Htx7RBuIf3iXCvmVZV4Rw5Ms5wb1JZr7O8umWFvid1fcx1XLkPhDvh
Pspwr5o2hMo95oOprn+5hHvI36zTx7FX+EMM96p9KzFFVwnhXvUt98ABwn3taUMI95i+1/Uv5oNk
nUq4jXBvepuyMRhiuKfctksI93XzgXCvmUa4pw+ENiv3tgMsRyl2mqYK5qb7anIK93ULBsI9AcI9
ru8xY9iE+7C0sT3sFrIepV431ulv6so9xXoZ2o9cw73Tk5gWL5NbtbJhp6oTgY4ciW8z9Boxi+8R
lzlurs91vo31CO1ZXFdiZVG57/7UCtkBcXH7ZVbufS2rDpV7c21Wu7HLpnKvzp8U35RHOSwT/iTD
X8SYD4iq9kLnW/W8Qo4iINyr2ygJ4R42rc9hmba3N8K94RtbF6rrnHRT/bzWP9RpVX8J9+XzlYRw
D5u2KgOanutCuCfQZ7g3XWFj50vxxjZtb502QhHu1WKOgiHc01TuMcsqOdwHd1XIEq8+16dczgat
29k4pJuEtHElTA4+QIzBVu4pKuH6/jY/ySa2vzFj86kriTptV+5Vr/3ur9p9fjCFLDv2vVwU+z7X
oXIfR+WeJNzN7KCk2zW7J+sd7v7HS+YZXLinXGFTP+ehhHvT/RFVy102ra9LCscEU+z6Ffp3dR8+
oae910l9XSfCPcNwN7M9kr4k6W2SnpX0qKSb3P30rvlGEe5NTzpK9drEVv9Nz/BblOK12d1e6HVL
Yt/3GCnCPeZ5xQZd6Dq1+H7dd590xRWzx+fPS4cPzx6nOqOWcM8z3N8iacPdb5z//xZJvrt6H0u4
h7S3Tn+7fG1SrLB1z7lKbOW+6rm0eQnkun6t298U7+U6/YvZVlK854R7t+GeYofqNZKeXvj/WUlv
TtAu0MjWjliz7UDf3Mx/p6y03ffFM4OX9XVIO5vRrY4vP7B54fF0OtGENXAwUl+KoKq9tlaJxaDf
3Nxe/gc+sD3ccPastHfv7HHd0EMX6l7fxeDf8vDDaY7MQQ6m2tycNm4l1bDMprsfnP+fYZmA+YY8
LFP3vEKnxbS5qr9149aLQzT79klPPTX798SJ5WPJR47sbK9uzPnAgbhT1pteZiPk98umMSzDsEyo
RyVdZ2bXSjon6SZJ70rQLrDSYhW7FcRblfbiN4PF+auq8N3fJhbb270j+vbbt+ffXUFvtb/sG8my
qruuT0CslIdC/qm2D4X8oyXzULmPoHJPfWOQut/Hvrcx2q5cm7ax6pvLkCr30MM461C5Jwr3oAUR
7qMI9xRSb+ip+5RruKdYH1KfdJViWur5CPfECPdywr3tM0P7XHbIslKfqJOijS7Xh7bCPeVtFwl3
wr2mv9VtxARMSeHetj6XXSXmgmDr/N2isYZ7m8VUVycXrvO8CPeWwz3kTU8RNiWHe4pqOpcLmDXV
xroSMo1wTzMf4R4h13Df2ce01X9V2+ssawjhXvc8x4ZwJ9yr2u863Ad3yV8gN3UnZOX6rWPZGbBI
I/UJf7E6rdw3NmbLir1KHZU7lfsYdFm5x7Qx9Mo9xRUzY7aVUQzLXDyNcCfcsYVw725YJrQfQwx3
hmUAjFouwyipUblH9CMGlTvq1A0VbP1u2bTFuz6VWrnHHk7a9JLPQ6/cCfeIfsRoI9xTn25OuA/X
qvcy5C5VuYZ7jNTrMuFet6AV4R5zVb0xhnvMWXyEe/n6DJ+6NuqmxZwBHIpwzyjcm3xNrJtWWri3
uxIR7kMV+16mvtdq7LSY+ULbGGu4Z79DtdSdHamtuvkFxzVjmRTH4g/xOP8xGFTlXjeNyj1s2TGV
Wuyy0J2215suv9VRuY9wWKZu2tDDPeSuPF1uVIT7sBDu1W10Ee7t7gsj3Acd7iHPk3BHFcK9uo0u
K/c2MotwJ9yDlx0yH+Gev1VDbG0d573Oe5z6ksehCHfCvbLNsYR76sPRCPc8tBmQbb/HhDvhnn24
p9h52Vc1TeU+bIR7+nBPfQtCwn3A4R6zrNj5QpdVp82bW6NbhHu7lXvT+er+jnAn3KOW1TbCPQ+E
O+Ge/UlM6+BEHQCYKbJyD22jbt6cKveURz20jcq9P6nvQ9vlPqi648FjULkT7hfNG3IBsxRCQzv1
it6GUm5ujZ3aLHbW2babtk+4t2wo4d60jVChbbZdJQFVCPdhh3tRY+7YRogjBfZjpdXlhRCp3Bv2
I1bblTvQ1OK61+Z4/jrTYtpvo3KPOfkvtvpnWIZwB5Jqc90beriH/E3stFThvmfdP9jZIXunmT1h
Zi+Z2f4mbQEA0mlUuZvZ6yW9LOkvJX3Q3R+rmZfKPbBNdoaiL12te1Tu4cvqdVjGzI5L+h3CPRzD
LRgzwj18WRwtE6iuMuly2dyODECbVoa7mR2TdNXiryS5pFvd/eg6C9tcOO5nMplo0kOiLQap2XbY
dr1sAFhuqs3NaeNWRj0s08bQDoDV2hiWafPeBEMclkkZ7h9098/VzEO4A5DU/ph7ivaGHu6NxtzN
7LCkP5d0paRPmNkJd7+xSZsAMEbLzgZuMow76pOYqNyBfgylcl/3iqwpKveLp3GGKuEODMRQwn3d
v8sp3BudoQoAyBPhDgAFItwBoECEOwAUiHAHgAIN+miZ48dXX8GOo2WA/HC0zDo3DB/5oZDVy40L
dy67C7QnNuhSLKvNv2vjQ4twr1xu88odQFpdbmvrLKvphwzhftG0tOEe/nWHcAf6kGu4t7kswj3g
96umhfeJcAf6kEvgdrkswj3g96umhfeJcAf6kEvgtr2sFPsSCPeoPhHuQB/GEu5p2uTaMgCAOcId
AApEuANAgbIbc489jDGuT4y5A30Y+jh4l8sqcodqfXuEOzBUQw/cLpc1uHCPOUQo9SnKhDvQj6EH
bpfLGly454BwB/ox9MDtclkcCgkAuGD0lXvIZYMBpDX0arrLZcVW7pek7cbwEOIAcrO4fzHW6Cv3
zLoEjMLQq+kul8UO1QiEO9CPtre9vm62Q7hngnAH+lHqtpdTuHO0DAAUiMo9ry4BxRrDfYlzqtwb
hbuZ3SbpHZK+Lek/Jb3H3b9ZMS/hDqBoOYV702GZByS9wd3fJOmMpA81bA8AkECjcHf3B9395fl/
H5G0t3mXAABNpdyh+l5Jn0rYHgAg0sozVM3smKSrFn8lySXd6u5H5/PcKulFd7+7lV4CANayMtzd
/Ya66WZ2s6RfkPRzq9ra3Ny88HgymWhSyi5yAEhkOp1q2vTaA2p+tMxBSX8i6Wfd/b9XzMvRMgCK
ltPRMk3D/YykSyVtBfsj7v4bFfMS7gCKVky4r7Ugwh1A4XIKdy4/AAAFItwBoECEOwAUiHAHgAIR
7gBQIMIdAArEoZB5dQnAwLR9nXqOc49AuAPIHce5AwAuINwBoECEOwAUiHAHgAIR7gBQIMIdAApE
uANAgUZ3nHvbJxwAQEqcxAQABeIkJgDABYQ7ABSIcAeAAhHuAFAgwh0ACkS4A0CBCHcAKBDhDgAF
ItwBoECEOwAUiHAHgAI1Cncz+30z+7yZPW5mnzazq1N1DAAQr2nlfpu7/5i7Xy/pHyRtJOgTAky3
Lm2Jxngt0+L1zEOjcHf3Fxb++12SXm7WHYRiA0qH1zItXs88XNK0ATP7A0m/Ium8pAONewQAaGxl
5W5mx8zs5MLPqfm/75Akd/+wu79W0t9J+s22OwwAWC3ZzTrM7DWSPunub6yYzp06ACBCzM06Gg3L
mNl17v7l+X8PS/pi1bwxnQMAxGlUuZvZvZJep9mO1K9K+nV3P5eobwCASJ3dQxUA0J3WzlA1s3ea
2RNm9pKZ7a+Z76CZnTazL5nZ77bVn6Ezs+81swfM7Ekz+0czu7xivpfM7LH5iWX3dd3PnK1a18zs
UjO7x8zOmNm/mdlr++jnUAS8nu82s6/P18fHzOy9ffRzCMzsDjN73sxO1szzZ/N184SZvWlVm21e
fuCUpF+U9HDVDGa2R9JHJP28pDdIepeZ/VCLfRqyWyQ96O6vl/SQpA9VzPe/7r7f3a9398PddS9v
gevar0r6hrv/oKTbJd3WbS+HY41t9575+rjf3f+m004Oy0c1ey2XMrMbJf3AfN38NUl/sarB1sLd
3Z909zOS6nakvlnSGXf/qru/KOkeSYfa6tPAHZJ05/zxnZrtwF6GHdfLhaxri6/xvZLe1mH/hiZ0
22V9DODu/yzpf2pmOSTprvm8/y7pcjO7qq7Nvi8cdo2kpxf+f3b+O1zsVe7+vCS5+3OSXlUx32Vm
9lkz+1cz44NyW8i6dmEed39J0nkze2U33Ruc0G33l+bDCB8zs73ddK1Iu1/vZ7QiK5seCnlM0uKn
h0lySbe6+9EmbY9Rzev54SWzV+0Jv9bdz5nZ90t6yMxOuvtXEnd1LKg6m7lf0t3u/qKZvU+zb0V8
G+pIo3B39xsaLv8ZSYs7rfbOfzdKda/nfGfLVe7+/Pzqm1+vaOPc/N+vmNlU0vWSCPewde2spNdI
etbMXiHpe9z9Gx31b2hWvp7uvjjM8NdiH0YTz2i2bm5ZmZVdDctUVUCPSrrOzK41s0sl3aTZpz0u
dr+km+eP3y3p47tnMLMr5q+jzOxKST8l6QtddTBzIevaUc1eW0n6Zc12XGO5la/nrkuAHxLr4iqm
6qy8X7NreMnM3iLp/NYwbSV3b+VHsx1+T0v6P0nnJH1q/vtXS/rEwnwHJT0p6YykW9rqz9B/JL1S
0oPz1+oBSVfMf//jkv5q/vgnJZ2U9Likz0u6ue9+5/SzbF2TdETS2+ePL5P0sfn0RyTt67vPOf8E
vJ5/KOmJ+fr4T5Je13efc/2RdLekZyV9W9LXJL1Hs6Ni3rcwz0ckfXm+be9f1SYnMQFAgfo+WgYA
0ALCHQAKRLgDQIEIdwAoEOEOAAUi3AGgQIQ7ABSIcAeAAv0/L+fz7Nq3eJAAAAAASUVORK5CYII=
"
>
</div>
</div>
</div>
</div>
</div>
<div class="cell border-box-sizing text_cell rendered">
<div class="prompt input_prompt">
</div>
<div class="inner_cell">
<div class="text_cell_render border-box-sizing rendered_html">
<h1 id="the-fitter">The fitter</h1>
<p>Here's a function to do the fitting. The actual calculation I defer to numpy's <code>lstsq</code> function, which uses the <a href="https://en.wikipedia.org/wiki/Singular_value_decomposition">singular value decomposition</a> under the hood and is supposed to be somewhat robust to degeneracies in the problem. (Not as robust as one could wish, unfortunately.) Because several of the error estimates are most easily calculated from some of these intermediate results, I calculate them here, though I will explain them below. Because we have quite a few values to return, I use a named tuple to return them.</p>
<p>(A brief aside: because scipy has only recently adopted the use of named return values, they don't want to break existing code that simply assigns the return values, as in <code>x, chi2, rk, s = lstsq(...)</code>, and named tuples satisfy this criterion. But I believe that in some cases they have objects that behave a little more like python's function call signatures: the named tuple is extended by more arguments accessible only by name. This would be nice here, but I don't know where to find the machinery to do it easily. Even the named tuple machinery is kind of cumbersome. I could write some, it wouldn't be long, but... I'll go down that rabbit hole another time.)</p>
<p>The basic fitting procedure can be carried out by matrix algebra. If we first leave aside the possibility of uncertainties, we get <span class="math display">\[ x = (A^TA)^{-1} A^T b. \]</span> Of course, numerically, matrix inversion is a troublesome operation, particularly when you have nearly-degenerate problems; the SVD provides a way to manage the treatment of such matrices, and the degree of near-degeneracy is summarized in the rank (<code>rk</code>) and singular-value (<code>s</code>) returns from <code>numpy.linalg.lstsq</code>. In such a numerically-messy case, though, the covariance matrix is going to become completely horrible, so we're going to ignore that when we're looking at error estimates. But in any case, we can work with the "<a href="https://en.wikipedia.org/wiki/Moore%E2%80%93Penrose_pseudoinverse">pseudo-inverse</a>", which provides some of the same resistance to numerical headaches (and is also based on the SVD). Anyway, incorporating uncertainties simply amounts to scaling the rows of <span class="math inline">\(A\)</span> and the elements of <span class="math inline">\(b\)</span> by dividing by the claimed uncertainties.</p>
</div>
</div>
</div>
<div class="cell border-box-sizing code_cell rendered">
<div class="input">
<div class="prompt input_prompt">In [3]:</div>
<div class="inner_cell">
<div class="input_area">
<div class="highlight"><pre><span class="n">linear_least_squares_result</span> <span class="o">=</span> <span class="n">namedtuple</span><span class="p">(</span><span class="s2">"linear_least_squares_result"</span><span class="p">,</span>
<span class="p">[</span><span class="s2">"names"</span><span class="p">,</span><span class="s2">"x"</span><span class="p">,</span><span class="s2">"chi2"</span><span class="p">,</span><span class="s2">"rk"</span><span class="p">,</span><span class="s2">"s"</span><span class="p">,</span><span class="s2">"res"</span><span class="p">,</span><span class="s2">"n"</span><span class="p">,</span><span class="s2">"cov"</span><span class="p">,</span><span class="s2">"cov_scaled"</span><span class="p">,</span><span class="s2">"leverage"</span><span class="p">,</span><span class="s2">"cov_robust"</span><span class="p">])</span>
<span class="k">def</span> <span class="nf">linear_least_squares</span><span class="p">(</span><span class="n">Acols</span><span class="p">,</span> <span class="n">b</span><span class="p">,</span> <span class="n">uncerts</span><span class="p">):</span>
<span class="n">names</span> <span class="o">=</span> <span class="nb">sorted</span><span class="p">(</span><span class="n">Acols</span><span class="o">.</span><span class="n">keys</span><span class="p">())</span>
<span class="n">A</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="n">Acols</span><span class="p">[</span><span class="n">n</span><span class="p">]</span> <span class="k">for</span> <span class="n">n</span> <span class="ow">in</span> <span class="n">names</span><span class="p">])</span><span class="o">.</span><span class="n">T</span>
<span class="n">n</span> <span class="o">=</span> <span class="nb">len</span><span class="p">(</span><span class="n">b</span><span class="p">)</span>
<span class="n">Ascaled</span> <span class="o">=</span> <span class="n">A</span><span class="o">/</span><span class="n">uncerts</span><span class="p">[:,</span><span class="bp">None</span><span class="p">]</span>
<span class="n">x</span><span class="p">,</span> <span class="n">chi2</span><span class="p">,</span> <span class="n">rk</span><span class="p">,</span> <span class="n">s</span> <span class="o">=</span> <span class="n">numpy</span><span class="o">.</span><span class="n">linalg</span><span class="o">.</span><span class="n">lstsq</span><span class="p">(</span><span class="n">A</span><span class="o">/</span><span class="n">uncerts</span><span class="p">[:,</span><span class="bp">None</span><span class="p">],</span> <span class="n">b</span><span class="o">/</span><span class="n">uncerts</span><span class="p">)</span>
<span class="n">res</span> <span class="o">=</span> <span class="n">b</span><span class="o">-</span><span class="n">np</span><span class="o">.</span><span class="n">dot</span><span class="p">(</span><span class="n">A</span><span class="p">,</span><span class="n">x</span><span class="p">)</span>
<span class="n">cov</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">linalg</span><span class="o">.</span><span class="n">pinv</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">dot</span><span class="p">(</span><span class="n">Ascaled</span><span class="o">.</span><span class="n">T</span><span class="p">,</span><span class="n">Ascaled</span><span class="p">))</span>
<span class="n">cov_scaled</span> <span class="o">=</span> <span class="n">cov</span><span class="o">*</span><span class="n">chi2</span><span class="o">/</span><span class="n">n</span>
<span class="n">leverage</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">diag</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">dot</span><span class="p">(</span><span class="n">Ascaled</span><span class="p">,</span><span class="n">np</span><span class="o">.</span><span class="n">dot</span><span class="p">(</span><span class="n">cov</span><span class="p">,</span><span class="n">Ascaled</span><span class="o">.</span><span class="n">T</span><span class="p">)))</span>
<span class="c1"># HC3, Mackinnon and White estimator</span>
<span class="n">cov_robust</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">dot</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">dot</span><span class="p">(</span>
<span class="n">cov</span><span class="p">,</span>
<span class="n">np</span><span class="o">.</span><span class="n">dot</span><span class="p">(</span><span class="n">Ascaled</span><span class="o">.</span><span class="n">T</span><span class="p">,</span><span class="n">Ascaled</span><span class="o">*</span><span class="p">(</span><span class="n">res</span><span class="o">/</span><span class="p">(</span><span class="mi">1</span><span class="o">-</span><span class="n">leverage</span><span class="p">))[:,</span><span class="bp">None</span><span class="p">]</span><span class="o">**</span><span class="mi">2</span><span class="p">)),</span>
<span class="n">cov</span><span class="p">)</span>
<span class="k">return</span> <span class="n">linear_least_squares_result</span><span class="p">(</span><span class="n">names</span><span class="o">=</span><span class="n">names</span><span class="p">,</span>
<span class="n">x</span><span class="o">=</span><span class="n">x</span><span class="p">,</span> <span class="n">chi2</span><span class="o">=</span><span class="n">chi2</span><span class="p">,</span> <span class="n">rk</span><span class="o">=</span><span class="n">rk</span><span class="p">,</span> <span class="n">s</span><span class="o">=</span><span class="n">s</span><span class="p">,</span> <span class="n">n</span><span class="o">=</span><span class="n">n</span><span class="p">,</span>
<span class="n">res</span><span class="o">=</span><span class="n">res</span><span class="p">,</span> <span class="n">cov</span><span class="o">=</span><span class="n">cov</span><span class="p">,</span> <span class="n">cov_scaled</span><span class="o">=</span><span class="n">cov_scaled</span><span class="p">,</span>
<span class="n">leverage</span><span class="o">=</span><span class="n">leverage</span><span class="p">,</span> <span class="n">cov_robust</span><span class="o">=</span><span class="n">cov_robust</span><span class="p">)</span>
</pre></div>
</div>
</div>
</div>
</div>
<div class="cell border-box-sizing code_cell rendered">
<div class="input">
<div class="prompt input_prompt">In [4]:</div>
<div class="inner_cell">
<div class="input_area">
<div class="highlight"><pre><span class="n">Acols</span> <span class="o">=</span> <span class="p">{</span><span class="s2">"constant"</span><span class="p">:</span> <span class="n">np</span><span class="o">.</span><span class="n">ones</span><span class="p">(</span><span class="nb">len</span><span class="p">(</span><span class="n">xs</span><span class="p">)),</span>
<span class="s2">"linear"</span><span class="p">:</span> <span class="n">xs</span><span class="p">}</span>
<span class="n">res</span> <span class="o">=</span> <span class="n">linear_least_squares</span><span class="p">(</span><span class="n">Acols</span><span class="p">,</span> <span class="n">ys</span><span class="p">,</span> <span class="n">true_uncerts</span><span class="p">)</span>
<span class="k">print</span> <span class="n">res</span><span class="o">.</span><span class="n">x</span>
</pre></div>
</div>
</div>
</div>
<div class="output_wrapper">
<div class="output">
<div class="output_area"><div class="prompt"></div>
<div class="output_subarea output_stream output_stdout output_text">
<pre>
[-0.00994284 -0.12111 ]
</pre>
</div>
</div>
</div>
</div>
</div>
<div class="cell border-box-sizing text_cell rendered">
<div class="prompt input_prompt">
</div>
<div class="inner_cell">
<div class="text_cell_render border-box-sizing rendered_html">
<h2 id="standard-linear-least-squares">Standard linear least squares</h2>
<p>This just uses the standard formula for the covariance matrix of a linear least squares fit: <span class="math display">\[ (A^T A)^{-1}, \]</span> for an <span class="math inline">\(A\)</span> with appropriately scaled rows. So it should work as long as it's provided with the correct uncertainties for the input values.</p>
</div>
</div>
</div>
<div class="cell border-box-sizing code_cell rendered">
<div class="input">
<div class="prompt input_prompt">In [5]:</div>
<div class="inner_cell">
<div class="input_area">
<div class="highlight"><pre><span class="n">res</span> <span class="o">=</span> <span class="n">linear_least_squares</span><span class="p">(</span><span class="n">Acols</span><span class="p">,</span> <span class="n">ys</span><span class="p">,</span> <span class="n">true_uncerts</span><span class="p">)</span>
<span class="k">print</span> <span class="n">res</span><span class="o">.</span><span class="n">names</span>
<span class="k">print</span> <span class="n">res</span><span class="o">.</span><span class="n">cov</span>
<span class="n">true_cov</span> <span class="o">=</span> <span class="n">res</span><span class="o">.</span><span class="n">cov</span>
</pre></div>
</div>
</div>
</div>
<div class="output_wrapper">
<div class="output">
<div class="output_area"><div class="prompt"></div>
<div class="output_subarea output_stream output_stdout output_text">
<pre>
['constant', 'linear']
[[ 4.15196328e-05 6.43500890e-19]
[ 6.43500890e-19 1.00000000e-02]]
</pre>
</div>
</div>
</div>
</div>
</div>
<div class="cell border-box-sizing text_cell rendered">
<div class="prompt input_prompt">
</div>
<div class="inner_cell">
<div class="text_cell_render border-box-sizing rendered_html">
<p>Let's quickly check that, statistically. That is, we want to know, if we had the same true values but many different instances of inaccurate measurements, what would the distribution of fit parameters be? It would be multivariate normal, described by a covariance matrix. We want our fitting algorithm to return that matrix, or an estimate of it. We can test whether it's right by actually creating many fake data sets, fitting them, and taking the covariance matrix of the fitted parameters. This should converge (like <span class="math inline">\(\sqrt{n}\)</span>) towards the true covariance matrix.</p>
<p>A note about how I structure the code: because the convergence requires many fake data sets, the more the better, I use three IPython blocks: one to set up the problem, one to run a bunch of simulations, and one to compute the results. They are written so I can rerun the middle block as many times as I have patience for and watch the computed results converge. If this gets slow, there are <a href="https://pythonhosted.org/joblib/parallel.html">tricks</a> to <a href="http://lighthouseinthesky.blogspot.nl/2014/10/parallel-ipython-notebooks.html">parallelize</a> the middle block, possibly even <a href="https://docs.python.org/3/library/concurrent.futures.html">running it in the background</a> while I evaluate partial results.</p>
</div>
</div>
</div>
<div class="cell border-box-sizing code_cell rendered">
<div class="input">
<div class="prompt input_prompt">In [6]:</div>
<div class="inner_cell">
<div class="input_area">
<div class="highlight"><pre><span class="k">def</span> <span class="nf">run_fake</span><span class="p">():</span>
<span class="n">ys</span> <span class="o">=</span> <span class="n">make_fake_data</span><span class="p">(</span><span class="n">xs</span><span class="p">,</span> <span class="n">true_uncerts</span><span class="p">)</span>
<span class="n">res</span> <span class="o">=</span> <span class="n">linear_least_squares</span><span class="p">(</span><span class="n">Acols</span><span class="p">,</span> <span class="n">ys</span><span class="p">,</span> <span class="n">true_uncerts</span><span class="p">)</span>
<span class="k">return</span> <span class="n">res</span><span class="o">.</span><span class="n">x</span>
<span class="n">fake_xs</span> <span class="o">=</span> <span class="p">[]</span>
</pre></div>
</div>
</div>
</div>
</div>
<div class="cell border-box-sizing code_cell rendered">
<div class="input">
<div class="prompt input_prompt">In [7]:</div>
<div class="inner_cell">
<div class="input_area">
<div class="highlight"><pre><span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="mi">100000</span><span class="p">):</span>
<span class="n">fake_xs</span><span class="o">.</span><span class="n">append</span><span class="p">(</span><span class="n">run_fake</span><span class="p">())</span>
</pre></div>
</div>
</div>
</div>
</div>
<div class="cell border-box-sizing code_cell rendered">
<div class="input">
<div class="prompt input_prompt">In [8]:</div>
<div class="inner_cell">
<div class="input_area">
<div class="highlight"><pre><span class="k">print</span> <span class="nb">len</span><span class="p">(</span><span class="n">fake_xs</span><span class="p">)</span>
<span class="n">fx</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="n">fake_xs</span><span class="p">)</span>
<span class="k">print</span> <span class="s2">"simulation:"</span>
<span class="k">print</span> <span class="n">np</span><span class="o">.</span><span class="n">cov</span><span class="p">(</span><span class="n">fx</span><span class="o">.</span><span class="n">T</span><span class="p">)</span>
<span class="k">print</span> <span class="s2">"claimed by linear_least_squares:"</span>
<span class="k">print</span> <span class="n">true_cov</span>
</pre></div>
</div>
</div>
</div>
<div class="output_wrapper">
<div class="output">
<div class="output_area"><div class="prompt"></div>
<div class="output_subarea output_stream output_stdout output_text">
<pre>
100000
simulation:
[[ 4.11844934e-05 4.03486016e-06]
[ 4.03486016e-06 1.00363575e-02]]
claimed by linear_least_squares:
[[ 4.15196328e-05 6.43500890e-19]
[ 6.43500890e-19 1.00000000e-02]]
</pre>
</div>
</div>
</div>
</div>
</div>
<div class="cell border-box-sizing text_cell rendered">
<div class="prompt input_prompt">
</div>
<div class="inner_cell">
<div class="text_cell_render border-box-sizing rendered_html">
<p>Well, they don't agree perfectly. The off-diagonal entries, which analytically should be exactly zero, are bad. But remember that error estimates converge like root n, and n here is not huge. So it looks like the errors here are okay.</p>
</div>
</div>
</div>
<div class="cell border-box-sizing text_cell rendered">
<div class="prompt input_prompt">
</div>
<div class="inner_cell">
<div class="text_cell_render border-box-sizing rendered_html">
<h2 id="dealing-with-wrong-error-estimates">Dealing with wrong error estimates</h2>
<p>Of course, we just saw the best-case scenario: the errors were normal and we knew exactly their sizes. In general we don't always know our uncertainties this well. What can be done about this? That is, suppose we put in our best guess for the experimental uncertainties and do the fitting. Is there any way to produce more pessimistic error estimates that might be okay even if our estimated uncertainties are wrong?</p>
<p>There's one commonly-used method: compute the reduced <span class="math inline">\(\chi^2\)</span> of the fit and rescale the uncertainties until that reduced <span class="math inline">\(\chi^2\)</span> is one. This correctly compensates for the situation where all your uncertainties are under- or over-estimated by the same factor. Of course, it suffers from the severe drawback that it quietly absorbs errors that are due to unmodelled structure (for example if these data had a quadratic part) into inflated errors on model parameters. If your model is wrong - that is, missing some important effect - then inflating your error bars is not at all a satisfactory solution, and it leaves you without much in the way of quality-of-fit estimates.</p>
<p>Nevertheless, what if we are confident that our model correctly describes all features of the data <em>but</em> that our uncertainties are all wrong by differing amounts? This situation - uncertainties that are different from each other - is called "<a href="https://en.wikipedia.org/wiki/Heteroscedasticity">heteroskedasticity</a>", and there is a fairly substantial statistical literature about how to deal with it.</p>
<p>First I'm going to repeat the trick above with the fake data to determine the "true" distribution of fitted parameters. It's going to be different from (worse than) the distribution above, because the fitting procedure is using the wrong weights in the least-squares fit.</p>
</div>
</div>
</div>
<div class="cell border-box-sizing code_cell rendered">
<div class="input">
<div class="prompt input_prompt">In [9]:</div>
<div class="inner_cell">
<div class="input_area">
<div class="highlight"><pre><span class="k">def</span> <span class="nf">run_fake_het</span><span class="p">():</span>
<span class="n">ys</span> <span class="o">=</span> <span class="n">make_fake_data</span><span class="p">(</span><span class="n">xs</span><span class="p">,</span> <span class="n">true_uncerts</span><span class="p">)</span>
<span class="n">res</span> <span class="o">=</span> <span class="n">linear_least_squares</span><span class="p">(</span><span class="n">Acols</span><span class="p">,</span> <span class="n">ys</span><span class="p">,</span> <span class="n">claimed_uncerts</span><span class="p">)</span>
<span class="k">return</span> <span class="n">res</span><span class="o">.</span><span class="n">x</span>
<span class="n">fake_het_xs</span> <span class="o">=</span> <span class="p">[]</span>
</pre></div>
</div>
</div>
</div>
</div>
<div class="cell border-box-sizing code_cell rendered">
<div class="input">
<div class="prompt input_prompt">In [10]:</div>
<div class="inner_cell">
<div class="input_area">
<div class="highlight"><pre><span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="mi">100000</span><span class="p">):</span>
<span class="n">fake_het_xs</span><span class="o">.</span><span class="n">append</span><span class="p">(</span><span class="n">run_fake_het</span><span class="p">())</span>
</pre></div>
</div>
</div>
</div>
</div>
<div class="cell border-box-sizing code_cell rendered">
<div class="input">
<div class="prompt input_prompt">In [11]:</div>
<div class="inner_cell">
<div class="input_area">
<div class="highlight"><pre><span class="k">print</span> <span class="nb">len</span><span class="p">(</span><span class="n">fake_het_xs</span><span class="p">)</span>
<span class="n">fhx</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="n">fake_het_xs</span><span class="p">)</span>
<span class="n">cov_het</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">cov</span><span class="p">(</span><span class="n">fhx</span><span class="o">.</span><span class="n">T</span><span class="p">)</span>
<span class="k">print</span> <span class="s2">"Observed covariance when using wrong uncertainties:"</span>
<span class="k">print</span> <span class="n">cov_het</span>
<span class="k">print</span> <span class="s2">"Observed covariance when using the right uncertainties:"</span>
<span class="k">print</span> <span class="n">true_cov</span>
</pre></div>
</div>
</div>
</div>
<div class="output_wrapper">
<div class="output">
<div class="output_area"><div class="prompt"></div>
<div class="output_subarea output_stream output_stdout output_text">
<pre>
100000
Observed covariance when using wrong uncertainties:
[[ 3.38980301e-03 -5.14196812e-05]
[ -5.14196812e-05 1.82065550e-02]]
Observed covariance when using the right uncertainties:
[[ 4.15196328e-05 6.43500890e-19]
[ 6.43500890e-19 1.00000000e-02]]
</pre>
</div>
</div>
</div>
</div>
</div>
<div class="cell border-box-sizing text_cell rendered">
<div class="prompt input_prompt">
</div>
<div class="inner_cell">
<div class="text_cell_render border-box-sizing rendered_html">
<p>Because we've weighted the data points uniformly, rather than based on their true uncertainty, we get an estimator that does a worse job than it could. I don't show it here, but it'll still converge to the right value eventually, if you keep adding data points, but it'll do so much more slowly than one using the right weights. So it pays to get your uncertainties as close to right as possible even if you're not going to trust them for error estimation.</p>
<p>Of course, if you make your covariance estimate with the wrong uncertainties, the covariance comes out badly wrong, in fact much too small in this case, and not even scaling so that the reduced <span class="math inline">\(\chi^2\)</span> is one solves the problem:</p>
</div>
</div>
</div>
<div class="cell border-box-sizing code_cell rendered">
<div class="input">
<div class="prompt input_prompt">In [12]:</div>
<div class="inner_cell">
<div class="input_area">
<div class="highlight"><pre><span class="n">res</span> <span class="o">=</span> <span class="n">linear_least_squares</span><span class="p">(</span><span class="n">Acols</span><span class="p">,</span> <span class="n">ys</span><span class="p">,</span> <span class="n">claimed_uncerts</span><span class="p">)</span>
<span class="k">print</span> <span class="s2">"Claimed covariance with the wrong uncertainties:"</span>
<span class="k">print</span> <span class="n">res</span><span class="o">.</span><span class="n">cov</span>
<span class="k">print</span> <span class="s2">"Claimed covariance with the wrong uncertainties, scaled:"</span>
<span class="k">print</span> <span class="n">res</span><span class="o">.</span><span class="n">cov_scaled</span>
<span class="k">print</span> <span class="s2">"Observed covariance when using wrong uncertainties:"</span>
<span class="k">print</span> <span class="n">cov_het</span>
</pre></div>
</div>
</div>
</div>
<div class="output_wrapper">
<div class="output">
<div class="output_area"><div class="prompt"></div>
<div class="output_subarea output_stream output_stdout output_text">
<pre>
Claimed covariance with the wrong uncertainties:
[[ 1.22500000e-03 1.27950983e-19]
[ -3.76251901e-19 3.60222772e-03]]
Claimed covariance with the wrong uncertainties, scaled:
[[ 2.92846580e-03 3.05877615e-19]
[ -8.99461898e-19 8.61142914e-03]]
Observed covariance when using wrong uncertainties:
[[ 3.38980301e-03 -5.14196812e-05]
[ -5.14196812e-05 1.82065550e-02]]
</pre>
</div>
</div>
</div>
</div>
</div>
<div class="cell border-box-sizing text_cell rendered">
<div class="prompt input_prompt">
</div>
<div class="inner_cell">
<div class="text_cell_render border-box-sizing rendered_html">
<h3 id="bootstrap">Bootstrap</h3>
</div>
</div>
</div>
<div class="cell border-box-sizing text_cell rendered">
<div class="prompt input_prompt">
</div>
<div class="inner_cell">
<div class="text_cell_render border-box-sizing rendered_html">
<p>The most direct approach to dealing with heteroskedasticity - and actually a pretty broad range of other kinds of statistical challenge - is the bootstrap method. In this method, fake - even more fake than the ones above - data sets are created from the observed data itself. Then the fitting method is applied repeatedly to these fake data sets to produce a large collection of fitted parameters, and their distribution is described. This has the advantage over the method above that one need know nothing about the method used to generate the true data, in particular, one need not know the true uncertainties. Annoyingly, it does take as many simulations as you have patience for (as opposed to the lovely matrix calculations you can do in one shot), but it's at least simple to implement.</p>
</div>
</div>
</div>
<div class="cell border-box-sizing code_cell rendered">
<div class="input">
<div class="prompt input_prompt">In [13]:</div>
<div class="inner_cell">
<div class="input_area">
<div class="highlight"><pre><span class="k">def</span> <span class="nf">run_boot</span><span class="p">():</span>
<span class="c1"># draw from original data with replacement</span>
<span class="n">ix</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">random</span><span class="o">.</span><span class="n">randint</span><span class="p">(</span><span class="nb">len</span><span class="p">(</span><span class="n">xs</span><span class="p">),</span> <span class="n">size</span><span class="o">=</span><span class="nb">len</span><span class="p">(</span><span class="n">xs</span><span class="p">))</span>
<span class="n">xs_l</span> <span class="o">=</span> <span class="n">xs</span><span class="p">[</span><span class="n">ix</span><span class="p">]</span>
<span class="n">ys_l</span> <span class="o">=</span> <span class="n">ys</span><span class="p">[</span><span class="n">ix</span><span class="p">]</span>
<span class="n">us_l</span> <span class="o">=</span> <span class="n">claimed_uncerts</span><span class="p">[</span><span class="n">ix</span><span class="p">]</span>
<span class="n">res</span> <span class="o">=</span> <span class="n">linear_least_squares</span><span class="p">({</span><span class="s2">"constant"</span><span class="p">:</span><span class="n">np</span><span class="o">.</span><span class="n">ones_like</span><span class="p">(</span><span class="n">xs_l</span><span class="p">),</span> <span class="s2">"linear"</span><span class="p">:</span><span class="n">xs_l</span><span class="p">},</span> <span class="n">ys_l</span><span class="p">,</span> <span class="n">us_l</span><span class="p">)</span>
<span class="k">return</span> <span class="n">res</span><span class="o">.</span><span class="n">x</span>
<span class="n">boots</span> <span class="o">=</span> <span class="p">[]</span>
</pre></div>
</div>
</div>
</div>
</div>
<div class="cell border-box-sizing code_cell rendered">
<div class="input">
<div class="prompt input_prompt">In [14]:</div>
<div class="inner_cell">
<div class="input_area">
<div class="highlight"><pre><span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="mi">100000</span><span class="p">):</span>
<span class="n">boots</span><span class="o">.</span><span class="n">append</span><span class="p">(</span><span class="n">run_boot</span><span class="p">())</span>
</pre></div>
</div>
</div>
</div>
</div>
<div class="cell border-box-sizing code_cell rendered">
<div class="input">
<div class="prompt input_prompt">In [15]:</div>
<div class="inner_cell">
<div class="input_area">
<div class="highlight"><pre><span class="k">print</span> <span class="nb">len</span><span class="p">(</span><span class="n">boots</span><span class="p">)</span>
<span class="n">ba</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="n">boots</span><span class="p">)</span>
<span class="k">print</span> <span class="s2">"Bootstrap:"</span>
<span class="k">print</span> <span class="n">np</span><span class="o">.</span><span class="n">cov</span><span class="p">(</span><span class="n">ba</span><span class="o">.</span><span class="n">T</span><span class="p">)</span>
<span class="k">print</span> <span class="s1">'"True":'</span>
<span class="k">print</span> <span class="n">cov_het</span>
</pre></div>
</div>
</div>
</div>
<div class="output_wrapper">
<div class="output">
<div class="output_area"><div class="prompt"></div>
<div class="output_subarea output_stream output_stdout output_text">
<pre>
100000
Bootstrap:
[[ 0.00294946 -0.00054469]
[-0.00054469 0.01402717]]
"True":
[[ 3.38980301e-03 -5.14196812e-05]
[ -5.14196812e-05 1.82065550e-02]]
</pre>
</div>
</div>
</div>
</div>
</div>
<div class="cell border-box-sizing text_cell rendered">
<div class="prompt input_prompt">
</div>
<div class="inner_cell">
<div class="text_cell_render border-box-sizing rendered_html">
<p>The bootstrap is doing a surprisingly bad job, here. The bootstrap might do better if we had more data points - a thousand rather than a hundred. Most of these techniques work only "asymptotically", that is, they converge towards being right as you supply more data points. The bootstrap is confusing because there are two numbers that go into it - the number of bootstrap trials and the number of data points. A small number of bootstrap trials has an obvious negative effect on the result, but even as the number of bootstrap trials goes to infinity, the number of data points still provides a limit on how well it can work.</p>
<p>There are smarter ways to apply the bootstrap to a linear least-squares problem; for example, maybe you could generate your fake data sets by adding residuals to the model, but with your residuals chosen with replacement from the set of residuals. This avoids the problem that some data points have more impact on the fitted parameters. This, by the way, is what the "<a href="https://en.wikipedia.org/wiki/Leverage_(statistics)">leverage</a>" measures; or rather, it measures the impact of the data points on the fitted model. There are <a href="https://en.wikipedia.org/wiki/Cook%27s_distance">derived measures</a> to quantify the impact of data points on fitted parameters more directly. But in case you're curious, here's a plot of the leverage of these data points.</p>
<p>The leverage is computed from the projection matrix (so called because it takes a data set to the best-fit model evaluated at those points) <span class="math display">\[ H = A(A^TA)^{-1}A^T \]</span> by picking out the diagonal entries <span class="math display">\[ h_i = \left[H\right]_{ii}. \]</span></p>
</div>
</div>
</div>
<div class="cell border-box-sizing code_cell rendered">
<div class="input">
<div class="prompt input_prompt">In [16]:</div>
<div class="inner_cell">
<div class="input_area">
<div class="highlight"><pre><span class="n">res</span> <span class="o">=</span> <span class="n">linear_least_squares</span><span class="p">(</span><span class="n">Acols</span><span class="p">,</span> <span class="n">ys</span><span class="p">,</span> <span class="n">true_uncerts</span><span class="p">)</span>
<span class="n">plt</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span><span class="n">xs</span><span class="p">,</span> <span class="n">res</span><span class="o">.</span><span class="n">leverage</span><span class="p">,</span> <span class="s2">"."</span><span class="p">)</span>
<span class="n">plt</span><span class="o">.</span><span class="n">ylabel</span><span class="p">(</span><span class="s2">"leverage"</span><span class="p">)</span>
<span class="n">plt</span><span class="o">.</span><span class="n">ylim</span><span class="p">(</span><span class="n">ymin</span><span class="o">=</span><span class="mi">0</span><span class="p">)</span>
<span class="n">plt</span><span class="o">.</span><span class="n">title</span><span class="p">(</span><span class="s2">"Using correct weights"</span><span class="p">)</span>
<span class="n">plt</span><span class="o">.</span><span class="n">figure</span><span class="p">()</span>
<span class="n">res</span> <span class="o">=</span> <span class="n">linear_least_squares</span><span class="p">(</span><span class="n">Acols</span><span class="p">,</span> <span class="n">ys</span><span class="p">,</span> <span class="n">claimed_uncerts</span><span class="p">)</span>
<span class="n">plt</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span><span class="n">xs</span><span class="p">,</span> <span class="n">res</span><span class="o">.</span><span class="n">leverage</span><span class="p">,</span> <span class="s2">"."</span><span class="p">)</span>
<span class="n">plt</span><span class="o">.</span><span class="n">ylim</span><span class="p">(</span><span class="n">ymin</span><span class="o">=</span><span class="mi">0</span><span class="p">)</span>
<span class="n">plt</span><span class="o">.</span><span class="n">ylabel</span><span class="p">(</span><span class="s2">"leverage"</span><span class="p">)</span>
<span class="n">plt</span><span class="o">.</span><span class="n">title</span><span class="p">(</span><span class="s2">"Using bogus weights"</span><span class="p">)</span>
</pre></div>
</div>
</div>
</div>
<div class="output_wrapper">
<div class="output">
<div class="output_area"><div class="prompt output_prompt">Out[16]:</div>
<div class="output_text output_subarea output_pyout">
<pre>
<matplotlib.text.Text at 0x7f64c095e450>
</pre>
</div>
</div>
<div class="output_area"><div class="prompt"></div>
<div class="output_png output_subarea ">
<img src="data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAY4AAAEKCAYAAAAFJbKyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz
AAALEgAACxIB0t1+/AAAH3xJREFUeJzt3X20XHV97/H3B5JYKyAqohIgAYMiwQei5lJbYCoKUcHU
pxqtFR+u1dJcV63eEm7rzYnWW9GuVi2lPqVe6qqN4LISUDBSGRAJIRIQyQMJCjEkEdfFQIJUCfFz
/9j7hMlwHmafzD5n5vB5rTVr9sNv//Z39pnZ37N/v/0g20RERHTqgIkOICIi+ksSR0REVJLEERER
lSRxREREJUkcERFRSRJHRERUksQRPUnS7ZJOneg4Jpsq21XSXZJeXndM0X+SOKI2kn4j6di2aYsl
fXm0ZW2faPu6+qLrLZKukfSuutfTre0q6TRJW7oRU/SfJI6o03BXl/b9VaeSHvPbGWraJCYmwd8x
xubx9EWP8acRZ0pPk3S5pB2S7pN0bcu8vc0k5VHKVyVdLGmnpB9JmtNSdo6kNZIekHSJpGWSPjLC
et8jaV1Z1+2SXlROP778z39HuY6zW5b5kqSLJH1T0i6gMcy0aZL+TtJmSdvL+U9oqWe+pFvKWDdJ
OkPS3wCnABeWMX1miJj/r6QPlMNHlEdzf1qOP1vSfS1lzyrXsUPS9ZKeP8x2/a1ym/5C0lpJ/3OI
o4iTJP2wrGtZ+fl+G/gWcISkXWXMz5T0Ukmry8+2XdLfjfT3j/6VxBET6YPAFuBpwOHA/xqh7NnA
V4AnA5cD/wQgaSrwdeBfgKcC/w68brhKJL0J+N/A22wfArwWuE/SlLLeq4CnA+8H/k3ScS2LvwX4
qO2Dge8PM+0CYBbwgvJ9erk+JM0FLgY+aPvJwKnA3bb/GvgesND2IbbfP0To1wKNcvg04Mfl8pTv
15XrOAlYCryn3B6fA5aX26ndAHA0MBN4JfA2HnsU8SbgDOCY8jO9w/ZDwKuAbbYPLmP+GfBp4FPl
Z3s2cMkQ64xJIIkjJtJu4FnAMbb32P7+CGWvt/1tFzdX+zLFTgzgd4ADbV9Y1vEfwE0j1PNu4BO2
1wDY/ontLcDJwJNsX2D7EdvXAFdQJIZBl9m+sVzu18NMew/wAdsP2P4l8PGWOt4FLLX93bL8dtsb
R9tIpWuB3yuHTwU+AfxuOX5aOZ9y/Z+1/QMXvgz8uvx87d4EfMz2TtvbgMcc6QCftn2v7fspEuuL
RojxYWCWpKfZfsj2SH+H6GNJHFGnPUD7f7pTKRIGwCcp/nNeIelOSeeNUNfPWoYfAn6r7FN4FrC1
rexInbZHletsd8QQy22mOGIYqd690yQ9Hfht4Oay+ecXwJUUR1QjrXtUtn8C/LI8ojiFIqltk/Qc
9k0cM4APDq5f0g7gyPLztTsCuGeUz3dvy/BDwEEjhPlu4LnABkmrJL2mg48WfSiJI+r0U4pmkFbH
UOyQsf2g7Q/ZfjZFk9FfSPr9iuvYzr47dyh20MPZQtGM0m7bEMsdzb5JaajO4NZp/49i5zrb9lPL
16Fl081I6x6u7nbXAm8EptreTtE8dQ5wKHBryzo+1rL+p9g+yPZXh6hvO0VSGXR0BzEMG6/tH9t+
q+2nUxwRfU3SEyvUGX0iiSPq9FXgryVNV+EVwFnApQCSXiNpcEe6C3iE4iilE4Md7yuBPZL+TNKB
kuYDc0dY7ovAhwY718uO5aOAVcBDkv5S0hRJjTLWf+/0w5bNaF8APlUefVB+9jPKIkuBd0r6/XJ7
HCHpueW8e4FjH1vrPq4DFpbvAM1y/Ho/+nyELwDvK/tTkPQkSa+W9KQh6rsEOF/SoZKmA3/W6Wct
432apEMGJ0j6I0mHlaMPUCSX31SoM/pEEkfU6SPADcD1wC8o2vvfant9Of844OryjKTvA//Uco3B
aP+BG8D2buD1wH8HdgBvpWiL//WQC9lfAz4GfEXSTuA/gKeW9ZwNvJriyOFC4I9tbxohnqGmnQfc
Cdwo6X5gBfCcct2rgXcCn6LYsTZ59L/8TwNvUnF22aeG+czXUjQVDTZLXQ88sWUc2zdT9HNcWDaV
baQ4Khkq5o9QHFHdVcZ5Kftut2H/BrbvoEiqPymbxJ4JzAPWltv1H4A3t/QFxSSiuh/kJGkexQ/l
AIqOwQuGKfcGii/uS2yvkTQDWA9sKIvcaPvcWoONSUHSjcA/2754omPpJ5LeR7Gzr9pcGI8zU+qs
vOy8vBA4naINebWky2xvaCt3EMXpjze2VXGn7TlEjEDFLTTuoDhSeBvwfIrTamME5VHCsRTNfc+h
OD16qDOrIvZRd1PVXGCT7c1lU8AyYP4Q5T5K0YzRflg74gVkEaXnAj+kaKr6APAG2/eOvEgA0yiu
89gJXE3RbPfPExpR9IVajzgoznZpPcXvHto6LsvTC4+0faWkv2xbfqakmym+2B+2fX2t0UZfsv0F
ik7hqMD2TymOziIqqTtxjEiSgL9n3867waOM7cDRtneUZ8B8Q9IJth8c7zgjIuJRdSeOrex7bviR
7Hte/MHAbKBZJpFnApdJem15Ze/DAGVn+Y8p2mHXtK5AUm60FhExBrbH1B1Qdx/HaopbEMyQNA1Y
ACwfnFne6uBw28faPoaic/zsMlEcVnauo+LW3LOAnwy1Ett5dem1ePHiCY9hMr2yPbM9e/W1P2o9
4rC9R9JCinPEB0/HXS9pCbDa9hXti/BoU9WpwEckPUxxEdF7XdwvJyIiJlDtfRy2r6I466V12uJh
yr68ZfjrFHc9jYiIHpIrx2MfjUZjokOYVLI9uyvbszfUfuV43SS53z9DRMR4k4R7tHM8IiImmSSO
iIioJIkjIiIqSeKIiIhKkjgiIqKSJI6IiKgkiSMiIipJ4oiIiEqSOCK6bNcuWLmyeB9qPKLfJXFE
dNGuXXDKKXDqqcX7tm37jid5xGSQxBHRRbffDmvXwiOPwLp18M1v7ju+du1ERxix/5I4IrroxBNh
9myYOhVOOAFe85p9x2fPnugII/ZfbnIY0WW7dhVHFrNnw8EHP3Y8ohfsz00OkzgiIh6HevruuJLm
SdogaaOk80Yo9wZJv5E0p2Xa+ZI2SVov6Yy6Y42IiNHV+gTA8pnhFwKnA9uA1ZIus72hrdxBwPsp
njk+OO15wB8CzwOOBK6WdFwOLyIiJlbdRxxzgU22N9veDSwD5g9R7qPAx4Fft0ybDyyz/Yjtu4FN
ZX0RETGB6k4c04EtLeP3lNP2knQScKTtK0dZdmv7shERMf5qbaoajSQBfw+csz/1DAwM7B1uNBp5
LnFERJtms0mz2exKXbWeVSXpZGDA9rxyfBFg2xeU44cAdwIPAgKeCdwHvBY4g6Lwx8uyVwGLba9q
W0e6PSIiKurls6pWA7MkzZA0DVgALB+caXun7cNtH2v7GIrO8bNtrynLvVnSNEnHALOAm2qONyIi
RlFrU5XtPZIWAisoktRS2+slLQFW276ifRGKIw9sr5N0CbAO2A2cm0OLiIiJlwsAIyIeh3q5qSoi
IiaZJI6IiKgkiSMiIipJ4oiIiEqSOCIiopIkjoiIqCSJIyIiKkniiIiISpI4IiKikiSOiIioJIkj
IiIqSeKIiIhKkjgiIqKSJI6IiKgkiSMiIipJ4oiIiEpqTxyS5knaIGmjpPOGmP9eSbdJukXSdZKO
L6fPkPSQpDXl66K6Y42IiNHV+gRASQcAG4HTgW0UzyBfYHtDS5mDbD9YDp9N8YjYV0maAVxu+wWj
rCNPAIyIqKiXnwA4F9hke7Pt3cAyYH5rgcGkUToI+E3L+Jg+VERE1KfuxDEd2NIyfk85bR+SzpV0
J/Bx4P0ts2ZKulnSNZJ+r95QIyKiE1MmOgAA2xcBF0laAHwYeAewHTja9g5Jc4BvSDqh7QgFgIGB
gb3DjUaDRqMxHmFHRPSNZrNJs9nsSl1193GcDAzYnleOLwJs+4JhygvYYfvQIeZdA3zQ9pq26enj
iIioqJf7OFYDs8ozpKYBC4DlrQUkzWoZPYuiMx1Jh5Wd60g6FpgF/KTmeCMiYhS1NlXZ3iNpIbCC
Ikkttb1e0hJgte0rgIWSXgE8DOwAzikXPxX4iKSHKTrM32v7/jrjjYiI0dXaVDUe0lQVEVFdLzdV
RUTEJJPEERERlSRxREREJUkcERFRSRJHRERUksQRERGVJHFEREQlSRwREVFJEkdERFSSxBEREZUk
cURERCVJHBERUUkSR0REVJLEERERlSRxREREJUkcERFRSe2JQ9I8SRskbZR03hDz3yvpNkm3SLpO
0vEt886XtEnSekln1B1rRESMrtYnAJbPDN8InA5so3gG+QLbG1rKHGT7wXL4bOBc26+SdALwb8BL
gSOBq4Hj2h/3lycARkRU18tPAJwLbLK92fZuYBkwv7XAYNIoHUTxfHGA1wLLbD9i+25gU1lfRERM
oCk11z8d2NIyfg9D7PwlnQv8BTAVeHnLsitbim0tp0VExASqO3F0xPZFwEWSFgAfBt5RZfmBgYG9
w41Gg0aj0cXoIiL6X7PZpNlsdqWuuvs4TgYGbM8rxxcBtn3BMOUF7LB9aHtZSVcBi22valsmfRwR
ERX1ch/HamCWpBmSpgELgOWtBSTNahk9i6IznbLcAknTJB0DzAJuqjneiIgYRa1NVbb3SFoIrKBI
Ukttr5e0BFht+wpgoaRXAA8DO4BzymXXSboEWAfspjjbKocWERETrNamqvGQpqqIiOp6uakqIiIm
mSSOiIioJIkjIiIqSeKIiIhKkjgiIqKSJI6IiKgkiSMiIipJ4oiIiEqSOCIiopIkjoiIqCSJIyIi
KkniiIiISpI4IiKiko4Th6QnSnpuncFERETv6yhxSDobuBW4qhx/kaTlIy8VERGTUadHHAPAXOB+
ANu3Asd0sqCkeZI2SNoo6bwh5n9A0lpJt0r6jqSjWubtkbRG0i2SvtFhrBERUaNOnwC42/YDxSPB
9xr16UmSDgAuBE4HtgGrJV1me0NLsTXAi23/StL7gE9SPGIW4Je253QYY0REjINOjzjWSnorcKCk
4yT9I3BDB8vNBTbZ3mx7N7AMmN9awPa1tn9Vjt4ITG+ZPaanU0VERH06TRz/A5gN/Br4d2An8Ocd
LDcd2NIyfg/7JoZ27waubBl/gqSbJN0gaf5wC0VExPjpqKnK9kPAX5WvWkh6G/Bi4LSWyTNsb5d0
DPBdSbfZvquuGCIiYnQdJQ5Jl/PYPo0HgB8An2tpamq3FTi6ZfzIclp7/a8AzgdOLZu0ALC9vXy/
S1ITOAl4TOIYGBjYO9xoNGg0GqN9pIiIx5Vms0mz2exKXbJH7eNG0qeBp1M0UwG8maK5ysAhtv94
mOUOBO6g6BzfDtwEvMX2+pYyJwGXAmfa/nHL9EOBh2w/LOkw4PvA/LaOdSS5k88QERGPkoTtMfUj
d3pW1ctsv7Rl/HJJq22/VNLa4RayvUfSQmAFRX/KUtvrJS0BVtu+AvgE8CTgUhWnbW22/QfA84DP
SdpTLvu37UkjIiLGX6dHHOspjgh+Wo4fDXzb9vMk3WL7pJrjHCm2HHFERFQ0HkccHwSul/RjilNk
jwHOlfQk4OKxrDgiIvpTR0ccAJKeABxfjt4xQof4uMoRR0REdftzxFElcZwInAD81uA02/86lpV2
UxJHRER1tTdVSVoMNCgSx7eAVwHXAxOeOCIiYnx1euX4GylOqf2Z7XcCLwSeXFtUERHRszpNHP9l
+zfAI5IOAX4OHDXKMhERMQl1elbVD8oL8r4A3Aw8CKysLaqIiOhZo3aOlxflHWl7Szk+k+Jq8dtq
j64D6RyPiKiu9rOqJP3I9vPHsoK6JXFERFS3P4mj0z6ONZJeOnqxiIiY7Do94tgAzAI2A7+kuHrc
tl9Qb3ijyxFHRER143HLkTPHUnlEREw+HTVV2d5Mcfrty8vhhzpdNiIiJpdOm6oWAy8Bnmv7OZKO
AC61/bt1BziaNFVFRFQ3Hp3jrwNeS9G/ge1twMFjWWFERPS3ThPHw+W/9QYob6ceERGPQ50mjksk
fQ44VNJ7gKspriIflaR5kjZI2ijpvCHmf0DSWkm3SvqOpKNa5p1TLneHpLd3GGtERNSoym3VXwmc
QXEq7rdtf6eDZQ4ANlLcIHEbsBpY0PoIWEmnAats/0rS+4CG7QWSngL8AJhTrvNmYI7tB9rWkT6O
iIiKxuO26n8BfLWTZNFmLrCpPBMLScuA+cDexGH72pbyNwJ/VA6fCawYTBSSVgDzgK9WjCEiIrqo
06aqg4EVkr4naaGkZ3S43HRgS8v4PeW04bwbuHKYZbeOsmxERIyDTq/jWGJ7NvBnwLOAayVd3c1A
JL0NeDHwyW7WGxER3dXpleODfg78DLgPOLyD8luBo1vGjyyn7UPSK4DzgVNt725ZttG27DVDrWRg
YGDvcKPRoNFoDFUsIuJxq9ls0mw2u1JXpxcAngv8IfB04FLgEtvrOljuQOAOis7x7cBNwFtsr28p
c1JZ55m2f9wyvbVz/IBy+MW2729bRzrHIyIqGo97VR0F/LntW6tUbnuPpIXACoqd/1Lb6yUtAVbb
vgL4BPAk4NLy2R+bbf+B7R2SPkqRMAwsaU8aEREx/qqcjvt7wHG2vyTp6cBBtu+qNbrO4soRR0RE
RePxIKfcqyoiYhLJvaoiImLc5F5VERFRSe33qoqIiMml1ntVjYf0cUREVFd753gvS+KIiKiutus4
JO2i7NdonwXY9iFjWWlERPSvEROH7Zw5FRER++i0czwiIgJI4oiIiIqSOCIiopIkjoiIqCSJIyIi
KkniiIiISpI4IiKikiSOiIiopPbEIWmepA2SNko6b4j5p0i6WdJuSa9vm7dH0hpJt0j6Rt2xRkTE
6Dp9dOyYSDoAuJDimePbgNWSLrO9oaXYZuAc4ENDVPFL23PqjDEiIqqpNXEAc4FNtjcDSFoGzAf2
Jg7bPy3nDXdPrIiI6CF1N1VNB7a0jN9TTuvUEyTdJOkGSfO7G1pERIxF3Ucc+2uG7e2SjgG+K+k2
23e1FxoYGNg73Gg0aDQa4xdhREQfaDabNJvNrtRV6/M4JJ0MDNieV44vorgd+wVDlP0ScLntrw9T
15Dz8zyOiIjq9ud5HHU3Va0GZkmaIWkasABYPkL5vR9C0qHlMkg6DHgZsK7OYCMiYnS1Jg7be4CF
wApgLbDM9npJSySdBSDpJZK2AG8EPivpR+XizwN+IOkW4D+Bv207GysiIiZAHh0bEfE41MtNVRER
MckkcURERCVJHBERUUkSR0REVJLEERERlSRxREREJUkcERFRSRJHRERUksQRERGVJHFEREQlSRwR
EVFJEkdERFSSxBEREZUkcURERCVJHBERUUkSR0REVFJ74pA0T9IGSRslnTfE/FMk3Sxpt6TXt807
p1zuDklvrzvWiIgYXa1PAJR0ALAROB3YRvEM8gWtj4CVdDRwCPAhYLntr5fTnwL8AJhD8Szym4E5
th9oW0eeABgRUVEvPwFwLrDJ9mbbu4FlwPzWArZ/avt2oH3vfyawwvYDtu+neG75vJrjjYiIUdSd
OKYDW1rG7ymnjWXZrRWWjYiImkyZ6AC6YWBgYO9wo9Gg0WhMWCwREb2o2WzSbDa7UlfdfRwnAwO2
55XjiwDbvmCIsl8CLm/p41gANGy/rxz/LHCN7a+2LZc+joiIinq5j2M1MEvSDEnTgAXA8hHKt36I
bwOvlPTksqP8leW0iIiYQLUmDtt7gIUUHdtrgWW210taIuksAEkvkbQFeCPwWUk/KpfdAXyU4syq
VcCSspM8IiImUK1NVeMhTVUREdX1clNVRERMMkkcERFRSRJHRERUksQRERGVJHFEREQlSRwREVFJ
EkdERFSSxBEREZUkcURERCVJHBERUUkSR0REVJLEERERlSRxREREJUkcERFRSRJHRERUksQR0QW7
dsHKlcV7N8pF9LLaE4ekeZI2SNoo6bwh5k+TtEzSJkkrJR1dTp8h6SFJa8rXRXXHGjEWu3bBKafA
qacW78MlhU7LRfS6WhOHpAOAC4EzgdnAWyQd31bs3cAvbB8HfAr4RMu8O23PKV/n1hlrxFjdfjus
XQuPPALr1hXD+1MuotfVfcQxF9hke7Pt3cAyYH5bmfnAxeXw14DTW+aN6bGGEePpxBNh9myYOhVO
OKEY3p9yEb1uSs31Twe2tIzfQ5FMhixje4+k+yU9tZw3U9LNwE7gw7avrzneiMoOPhi+973iCGL2
7GJ8f8pF9Lq6E8dYDB5lbAeOtr1D0hzgG5JOsP1g+wIDAwN7hxuNBo1GYzzijNjr4IPh5JO7Vy6i
25rNJs1msyt1yXZXKhqyculkYMD2vHJ8EWDbF7SUubIss0rSgcB224cPUdc1wAdtr2mb7jo/Q0TE
ZCQJ22PqDqi7j2M1MKs8Q2oasABY3lbmcuCccvhNwHcBJB1Wdq4j6VhgFvCTmuONiIhR1NpUVfZZ
LARWUCSppbbXS1oCrLZ9BbAU+LKkTcB9FMkF4FTgI5IeBn4DvNf2/XXGGxERo6u1qWo8pKkqIqK6
Xm6qioiISSaJIyIiKkniiIiISpI4IiKikiSOiIioJIkjYgy6cXv03GI9+lUSR0RF3bg9em6xHv0s
iSOiom7cHj23WI9+lsQRUVE3bo+eW6xHP8uV4xEd2rWrOFI48cRifH9vj75r1751tNafW65H3fbn
yvEkjogRDO7MZ8yAV7/60R39977X3Z37YJ/HYP3f+hZs3pwkEvVJ4ujzzxC9o/2oYnBnPnMm3H13
0ScxdSpcd113n6uxcmXRUf7IIzBlyqPrSxKJuuReVfG41XpKa/vprcPNG67ctm37num0atWjHdh3
313szOvqk2jt82hNUmvXwmmn7Xv2VTc/c87mijGx3dcvwDt32ra9c6d9ww3Fe+vwSPM6LdeNOibr
uiYq3q1b7Re+0J4yxT7xxOI1ZUoxbbh5I5WbNat4B3vqVPvqq4t5U6c+Wnblykfj6LadO4v6B2Oa
OnX4mLr1mQfHJ/pv2Q/r6rd4RytX7P7HuN8d64K98gJG/NGM9QfV7Tom67omMt7WneqBB+67g/38
54eeN1K5wToHE8XgD6zOZDGcoZLIC19of+c79Xzmif5b9vq6+i3eTsr1dOIA5gEbgI3AeUPMnwYs
AzYBKymeMz447/xy+nrgjGHqH/FHM9YfVLfrmKzrmsh4W3f0gz+G1qODwR1u67yRyo3HUcVYtCav
nTuH/lxj/cztRzT99N17vHzP61tXjyYOij6UO4EZwFTgVuD4tjJ/ClxUDr8ZWFYOnwDcQvGUwpll
PRpiHSP+aMb6g+p2Hf2yrpkzr+mbeNt39O1HB+073E7Kdds111zT9To7/SydlBtqe/by93zmzGsm
5DfVD9umarleThwnA1e2jC9qP+oArgL+Wzl8IPDzocoCVw6Wa1t+1B/NWH5QddTRD+tatGhxX8Xb
6xYvXjzRIYyqV/6WnZRbtGjxhP2men3bVC3Xy4njDcDnW8bfBnymrcyPgCNaxjcBTwX+EXhry/Qv
Aq8fYh3D/BxiLPphR9dPsj27K9uze/YncfTi6bhjOq84IiLGR60XAEo6GRiwPa8cX0SR5S5oKXNl
WWaVpAOB7bYPby8r6Spgse1VbevI1X8REWPgMV4AOKXbgbRZDcySNAPYDiwA3tJW5nLgHGAV8Cbg
u+X05cC/SfoHYDowC7ipfQVj/eARETE2tSYO23skLQRWUJxhtdT2eklLgNW2rwCWAl+WtAm4jyK5
YHudpEuAdcBu4FzXeXgUEREd6ft7VUVExPjqxc7xEUl6o6TbJe2RNGeEcvMkbZC0UdJ54xljP5H0
FEkrJN0h6duSnjxMuT2S1ki6RdI3xjvOXjfa903SNEnLJG2StFLS0RMRZz/oYFueI+nn5fdxjaR3
TUSc/UDSUkn3SrpthDKfKb+Xt0p6USf19l3ioDh993XAtcMVkHQAcCFwJjAbeIuk48cnvL6zCLja
9nMp+pfOH6bcL23PsX2S7T8Yv/B6X4fft3cDv7B9HPAp4BPjG2V/qPDbXVZ+H+fY/pdxDbK/fIli
Ww5J0quAZ5ffy/cCn+2k0r5LHLbvsL2JkU/bnQtssr3Z9m6KW5rMH5cA+8984OJy+GJguKSQkxCG
18n3rXU7fw04fRzj6yed/nbzfeyA7euBHSMUmQ/8a1l2FfBkSc8Yrd6+Sxwdmg5saRm/p5wWj3W4
7XsBbP8MOHyYck+QdJOkGyQlCe+rk+/b3jK29wD3S3rq+ITXVzr97b6+bFq5RNKR4xPapNS+vbfS
wb6y7tNxx0TSd4DWrCfAwF/ZvnxioupfI2zPvx6i+HBnS8ywvV3SMcB3Jd1m+64uh/p4kv+Yx245
8BXbuyX9CcWRXI7gxlFPJg7br9zPKrYCrZ2PR5bTHpdG2p5lx9kzbN8r6ZnAz4epY3v5fpekJnAS
kMRR6OT7dg9wFLCtvND1ENu/GKf4+smo29J2a9PLF0l/0f7YSvG9HNTRvrLfm6qG+69t74WHkqZR
XBuyfPzC6ivLgXeUw+cAl7UXkHRouR2RdBjwMorra6LQyfdt8EJX2PdC19jXqNuy/Adn0HzyXRyN
GH5fuRx4O+y908f9g03XIxrrTa4m6kXRebsF+C+Kq9GvLKc/C7iipdw84A6KmyYumui4e/VFcUPJ
q8tttQI4tJz+YsobVAK/A9xGcZv7HwLvmOi4e+011PcNWAKcVQ4/AbiknH8jMHOiY+7VVwfb8v8A
t5ffx/8EnjPRMffqC/gKsA34NfBT4J0UZ0/9SUuZCykeW/FDYE4n9eYCwIiIqKTfm6oiImKcJXFE
REQlSRwREVFJEkdERFSSxBEREZUkcURERCVJHBERUUkSR0REVPL/AYVM1mDwLGpKAAAAAElFTkSu
QmCC
"
>
</div>
</div>
<div class="output_area"><div class="prompt"></div>
<div class="output_png output_subarea ">
<img src="data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAZQAAAEKCAYAAAA1qaOTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz
AAALEgAACxIB0t1+/AAAIABJREFUeJzt3X+8HHV97/HXOwkBxSBGRSSYQCWCJv4ASkitlaMoREsN
tiKJ15paqreita0/Gmy9l5PeR6vQ9haV1nrbaJHaRopaAfklwjEokRwFQUgCp0pCfoFICAkqmoTP
/WNmk8lm9pzdc3Z2Z2ffz8djH5md+e7Mdzc738/5/lxFBGZmZhM1qdsZMDOzanBAMTOztnBAMTOz
tnBAMTOztnBAMTOztnBAMTOztnBAsZ4m6R5Jry7gvEsk3dru85aNpGsl/W6TaW+R9PtF58l615Ru
Z8BM0lPAcRHxo8y+C9N9oxZ2ETG3wKxVfpJWRLyxHeeRNAt4AJgSEU+145zWe1xDsTJoVHBXvkCv
EJH8f6nbGbHucUCxMhi1EJL0bElXS3pM0qOSvpk59oCk16bbF0r6oqTLJO2Q9ANJJ2XSniTpDkmP
S7pC0gpJfznKpSdJ+pSk7ZLW1K6Tnuv5kr6a5ud+SX+QOXZImodtku6V9GFJGzPHn5L0K5nnn6vl
Y7T3WveZDEr6ZLo9RdITki7KXP/nkg5Pn8+X9O30nHdKOi1znr3NWJImSfo7SY9I+qGk96Z5zZYT
x0j6Vvr5Xi9perq/ls/t6bFTJb1Q0lD6+f1Y0n+M8llbBTigWC/4ILAReDZwBPDno6T9LeDfgWcC
VwP/ACDpIODLwGeB6cB/AG8e47qnAiPpdQeBL9cKaeCLwIPAkcA5wF9LGkiPDQIzgWOA1wNvZ//a
1mg1r2bf6zeBWmA4BXgIqPUlvRJYFxHbJc0ArgH+MiKeBXwI+JKkZ+ec893AmcDLgJOAs3PyuhhY
AjwXODg9H5lrHxYRh0XE7cD/AW6IiMOBo4FPjfK+rQIcUKwX7AKeDxwbEXsi4tujpP1WRNwQySJ1
l5MUjgC/BkyOiEvTc3wFWD3GdR+OiE+m6a8A7gN+U9LR6fmWRsSuiLgL+BfgHenrzgH+KiJ2RMQW
4JN15x2tRtbse10FzJb0LJLCfDkwQ9LT0+e1GsP/AL4WETcARMQ3gO8CeX0n5wCfiIitEfE48PGc
NJ+LiB9GxC+AK4BXjPLedgGzJM2IiF9GxG2jvG+rAAcUK4M9wEF1+w4iKZAA/gb4IXCjpP+WtHSU
cz2U2f4ZcEjaZPN8YHNd2o2Mrj79BuCo9LEtIn5Wd2xGun0UsKmF62Q19V4j4kmSwDBAEkCGgNuA
V5HUXGoBZRbw1rT5bZukx4BfJ6lZ1TuqLq95+a7/fJ8xynv5MEkZszptfnznKGmtAhxQrAweJGke
yjqWpJAmIp6IiA9FxAuBNwEfkPSaFq+xlX0Ffs0LxnhNffqZwJb0MV3SoXXHagFoK0kTT/ZY1s+A
p2ee7y3cW3yvK4HXktQShtPnZ5I0ga1M02wEPh8R09PHsyJiWkT8Tc75xsr3aA5oxouIH0fEuyNi
BvCHwD9m+46sehxQrAy+CHxU0gwlXgecBfwngKTflPTCNO1OYDdJraYZtSaYVcCetKN5sqSFwLwx
Xvs8SX+UdnqfA5xA0ny0iaQ28DFJB0t6GXAeSRMbJE1BH5F0eNqH8d66894JvC3tBF/Avr6QRu+1
0TDcb5I0s62JiN0ktZQ/AB6IiEfTNP8G/JakM9LrHSLpNElH5ZzvCuCPJR2V9hX92RifT9YjaT5r
eUfSW9L3D7A9Pe4hxRXmgGJl8JckBfS3gG0kbfdvi4i16fHZwE2SdgLfBv4hImp/gY81tDgAImIX
8NskBe5jwNtIOu1/Mcprv5Ne+yckHcy/ExHb02OLSWpRW4AvAf8rIm7JvJ/NJPMybiQJjNnr/AlJ
7eOx9DxfyRzLe6+5I71IPrNDSJu3ImIN8HP2NXeRBr+FJJ37j5DU+j7Evns/+/n9c5rfu4HvAV8D
dmfmlTT8rCPi58BfAd9Om9bmkdSUbpe0A/gv4P0Rsb7ROaz3qegf2Er/AruE5Au8PCIuqjs+Ffg8
cDLJjXtuRDyYOT4TuBe4MCL+bzPnNGuGpO8An46Iywq+zh+SfK9bbabrqvQ++3REHNvtvFhvKLSG
knaGXkrSrjsHWCzphLpk55F0cM4mCRIX1x3/O+DaFs9pdgBJr5b0vLTJawnwUuD6Aq5zpKRXps13
x5MMBf5yu6/Tbmlz2BvSz2cGcCE9kG8rj6KbvOYBIxGxIW1yWEFS/c5aCNT+QrwSOL12IG3n/hFJ
DaWVc5rlOR64i6Sp6U9JmrAeLuA6U4HPADuAm0iatD5dwHXaTcAykmbH75G2DHQ1R9ZTil7Lawb7
Dz3cxIEdoXvTRMSedFbtdJI25z8jmRj24RbPaXaAiPhnkn6Coq/zIEntp6ek/SC+l2zcytgpXxuV
Mwj8fd1YfzMzK6miayib2X8s+9EcOFlsE8l8gC2SJpMs3bBN0qnA70i6GHgWyZDPJ4E7mjgnAJK8
uKCZ2ThERMsLfRZdQxkGjpM0Kx3NtQi4qi7N1SRrA0Gy9MPNABHx6oj4lYj4FZLO+r+OiH9s8px7
7dgRRPgx0ceFF17Y9TxU6eHP059nWR87doz/7/BCA0pE7AHeRzK2/V5gRUSslbRM0llpsuXAcySN
kIzPv2A852yU/jd+A3bunPh7MTOrup07kzJzvAr/ga2IuJ5kdE1234WZ7V8Abx3jHMvGOmcja9bA
vffC/PlNZ9nMrC/dc09SXo5XGTvl2+olL4E5c7qdi943MDDQ7SxUij/P9vLn2R5z506svCx8pnw3
SYodO4Jp07qdEzOz3rBzJxx2mIhxdMpXPqBU+f2ZmRVBGl9AqXyTFyQRd9Uqd86bmdVrZ/lY+YBS
G7Xw6ld7xJeZWVa7y8fKB5TaqIXdu/eN+DIzs/aXj5UPKLVRCwcd5BFfZmZZ7S4f+6JTfufOJPLO
mYNHfJmZZeSVj+PtlO+LgGJmZs3zKC8zM+uqvgsoHkJsZv2uqHKwrwKKhxCbWb8rshzsq4DiIcRm
1u+KLAf7KqB4CLGZ9bsiy8G+G+XlIcRm1u/GKgc9bDiHhw2bmbXOw4bHyaO+zKzqOlXO9XVA8agv
M6u6TpZzfR1QPOrLzKquk+Vc4QFF0gJJ6yTdL2lpzvGpklZIGpG0StLMdP8pku7MPM7OvGa9pLvS
/avHmzeP+jKzqutkOVdop7ykScD9wOnAFmAYWBQR6zJp3gO8NCLOl3Qu8OaIWCTpEOCXEfGUpCOB
u4Dnp89/BJwcEY+Ncf0xO+U96svMqq7Vcq6snfLzgJGI2BARu4AVwMK6NAuBy9LtK0mCDxHxZEQ8
le5/GvBU5jWiTXmfNg3mz3cwMbPq6lQ5V3RAmQFszDzflO7LTRMRe4DtkqYDSJon6R6S2skfZgJM
ADdIGpb0rnZl1iO+zKwqulGelbFTfm81KyJWR8Rc4BTgzyVNTQ/9ekT8KvBG4L2SXjXRi3rEl5lV
RbfKsykFn38zMDPz/Oh0X9Ym4AXAFkmTgcMiYls2QUTcJ+kJYC5wR0RsTfc/IukrJE1r38rLwODg
4N7tgYEBBgYGcjOaNxJi/vym36eZWWm0Wp4NDQ0xNDQ04esW3Sk/GbiPpF9kK7AaWBwRazNpzgfm
pp3yi4Cz0075Y4CNEbFH0izg28DLgCeBSRHxhKRDgRuBZRFxY871m54pX4voa9YkIyFuvdX9KmbW
myZanpV26RVJC4BPkDSvLY+Ij0taBgxHxDWSDgYuB04EHiUZBbZe0tuBC4BfknTIL4uIqyUdC3yF
pB9lCvCFiPh4g2u3tPSKR3yZWVVMpDwrbUDpJq/lZWbWurIOG+5pHvVlZr2k22WWA0oDHvVlZr2k
DGWWA0oDXufLzHpJGcosB5QGvM6XmfWSMpRZ7pQfhUd9mVkvaVeZ5VFeOTzKy8ysdR7lVbBuj54w
M6tXtnLJAaUJZRg9YWaWVcZyyQGlCWUYPWFmllXGcskBpQllGD1hZpZVxnLJnfJN8ogvMyubosol
j/LK4VFeZmat8yivDirbyAoz6x9lLn8cUFpUxpEVZtYfyl7+OKC0qIwjK8ysP5S9/HFAaVEZR1aY
WX8oe/njTvlx8IgvM+uWTpQ/HuWVw6O8zMxa51FeXVTmURdm1tt6qXwpPKBIWiBpnaT7JS3NOT5V
0gpJI5JWSZqZ7j9F0p2Zx9nNnrOTyj7qwsx6V6+VL4UGFEmTgEuBM4E5wGJJJ9QlOw/YFhGzgUuA
i9P9PwBOjogTgTcAn5E0qclzdkzZR12YWe/qtfKl6BrKPGAkIjZExC5gBbCwLs1C4LJ0+0rgdICI
eDIinkr3Pw2obTdzzo4p+6gLM+tdvVa+TCn4/DOAjZnnm0gCQm6aiNgjabuk6RGxTdI84LPATOB3
I+IpSc2cs2OmTYNbb/WoLzNrv14rX4oOKOOxd2RBRKwG5ko6Hvi8pOtaPdng4ODe7YGBAQYGBtqQ
xf1Nmwbz5yfbO3cm1dS5c8v/n29m5VRfjtTKl6IMDQ0xNDQ04fMUOmxY0nxgMCIWpM8vACIiLsqk
uS5Nc7ukycDWiDgi51zfAD4MTB3rnJnXdHTYcK0DrfbXxK23OqiYWWvKUI6UddjwMHCcpFmSpgKL
gKvq0lwNLEm3zwFuBpB0TBpgkDQLOB5Y3+Q5u6LXOtDMrHx6uRwpNKBExB7gfcCNwL3AiohYK2mZ
pLPSZMuB50gaAf4EuCDd/yrgLkl3AF8C3hMR2xqds8j30axe60Azs/Lp5XLEM+XbzMuymNlEdbsc
8dIrObz0iplZ68rah9L3emnZBDPrniqUFQ4oBeq1ZRPMrDuqUlY4oBSol0drmFnnVKWscEApUC+P
1jCzzqlKWeFO+YJ1e7SGmfWGMpUVHuWVowwBJcvLsphZTZnLA4/yKrmqdLqZ2cRVtTxwQOmQqnS6
mdnEVbU8cEDpkKp0upnZxFW1PHAfSgeVqdPNzLqrzOWBO+VzlC2gZJW5Q87MitEr97075XtIVTvk
zKyxfrjvHVC6oKodcmbWWD/c9w4oXVDVDjkza6wf7nv3oXRJmTvkzKwYvXLfu1M+R5kDSr1e6awz
s9b04r3tTvke1g+ddWb9qN/ubQeUEuiHzjqzftRv93bhAUXSAknrJN0vaWnO8amSVkgakbRK0sx0
/+skfVfSXZKGJb0m85pb0nPeKekOSc8p+n0UqR8668z6Ub/d24X2oUiaBNwPnA5sAYaBRRGxLpPm
PcBLI+J8SecCb46IRZJeDjwcEQ9JmgPcEBFHp6+5BfhARNw5xvV7qg+lFzrrzKw1vXhvl7UPZR4w
EhEbImIXsAJYWJdmIXBZun0lSfAhIu6KiIfS7XuBQyQdlHldpZrrpk2D+fOTf6vw29Jm/ar+/s3e
21VXdKE8A9iYeb4p3ZebJiL2ANslTc8mkPQW4I40KNV8Nm3u+mj7s909/daJZ1Yl/X7/Tul2BnLs
V81Km7s+Brw+s/ttEbFV0qHAlyW9PSL+Le9kg4ODe7cHBgYYGBhoe4bbKa8Tb/78bufKzJrRq/fv
0NAQQ0NDEz5P0X0o84HBiFiQPr8AiIi4KJPmujTN7ZImA1sj4oj02NHAN4AlEfGdBtdYApwcEe/P
OdYzfSg1tb9w1qxJOvFuvbU/qspmVVCV+7esfSjDwHGSZkmaCiwCrqpLczWwJN0+B7gZQNLhwDXA
0mwwkTRZ0rPT7YOAs4B7Cn0XHTRtWvIlXLky+Rfcn2JWdrV+E9j//u3FYDIRhc+Ul7QA+ARJ8Foe
ER+XtAwYjohrJB0MXA6cCDxKMgpsvaS/AC4ARkiawQI4A/gZsJKkuW4ycBPJiK8D3kgv1lCyan/t
1EaI9OMX1KzsqnifeumVHL0eUFatSjr3du9OxrGvXNkb7bFm/aSK92lZm7xsAvptUpRZL/J9uo9r
KCXXi5OizPpN1e5TN3nlqEJAqdeLK5eaVU3V70M3efWBfp80ZVYGvg8bazqgSHqapOOLzIyNrt9W
LjUrI9+HjTUVUCT9FvB94Pr0+Ssk1c8nsYK588+s+3wfNtZUH4qk7wGvBYYi4sR03w8i4qUF529C
qtqHUuv8g2q345qVSbbfBKrVCV+v6D6UXRHxeN2+apXUPaK2cim4HdesU+r7TaB/VhBuRbMB5V5J
bwMmS5ot6VPAbQXmy8bgdlyzzvH91pxmA8ofAXOAXwD/AewA/qSoTNnY3I5r1jm+35rjeSg9rGqT
qczKrJ/ut0InNkq6mgP7TB4Hvgt8JiKebPXCnVD1gFKv6pOtzDqtX++pojvlfwQ8Afxz+tgB7ARe
lD63LvNkK7P28j3VumZ/sfGVEXFK5vnVkoYj4hRJ7p4qgV79pTizsvI91bpmayjPkDSz9iTdfkb6
9Jdtz5W1zJ2GZu3le6p1zfahvBH4J+CHJD92dSxwPjAEvCsiLikwj+PWj30onvRoNn71fSb91BGf
Vfhqw+kvK56QPr2vrB3xWf0WUGqq+AtyZkXzfbNPJ1Ybng0cD7wceKukd7R6MesMT8Iya53vm4lr
dnHIC4FPpY/XABcDbyowXzYBbvs1a53vm4lrtobyFuB04KGIeCdJLeWZzbxQ0gJJ6yTdL2lpzvGp
klZIGpG0qtb5L+l1kr4r6S5Jw5Jek3nNSZLuTs9Zyv6bbpo2Lamur1yZ/AvJ71572KPZgXbuTO4P
2P++6dfmroloNqD8PCKeAnZLOgz4MfCCsV4kaRJwKXAmydItiyWdUJfsPGBbRMwGLiGp/QA8ApwV
ES8Hfg+4PPOaTwPnRcSLgBdJOrPJ99E3vIik2di86GN7NRtQvivpcJJJjN8D7gBWNfG6ecBIRGyI
iF3ACmBhXZqFwGXp9pUkNSEi4q6IeCjdvhc4RNJBko4EpkXEcPqazwNnN/k++o7bhc0a8/3RXmMG
FEkCPhYR2yPin4DXA0vSpq+xzAA2Zp5vSvflpomIPcB2SdPr8vAW4I40KM1IzzPaOS3ldmGzxnx/
tNeYM+UjIiRdC7w0fb6+4DztN1RN0hzgYySBrGWDg4N7twcGBhgYGJhA1npPrT8lO5a+X9cnMoMD
v//190c/GhoaYmhoaMLnaXZi42XApZlmpuZOLs0HBiNiQfr8ApIYdVEmzXVpmtslTQa2RsQR6bGj
gW+Q1Ii+k+47ErglIl6cPl8EnBYR78m5fl/OQxmNx9pbP/P3vzlFz0M5FVgl6Yfp6KofSLq7idcN
A8dJmiVpKrAIqP8t+quBJen2OcDNAGmfzTXA0lowAUj7VR6XNC9tjnsH8NUm30ffc5ux9TN//4vV
7OKQ4xpFFRF7JL0PuJEkeC2PiLWSlgHDEXENsBy4XNII8ChJ0AF4L/BC4H+n82ACOCMifpIe+1fg
EODaiLh+PPnrR7U24zVrkjbjmTOTIZNu/rIqqzVzzZq1//fffSbt1crSK68CZkfE5yQ9F3hGRDxQ
aO4myE1e+WrrE82cCW98o6v/Vm31zVzXXgsPPtjffSZjKbTJK60hLAU+ku46CPi3Vi9m5VCbo7Jh
g6v/Vn31zVwPPui5JkVptg/lzSRLrfwUICK2AP7v6HEeMmn9wN/zzmk2oPwybTsKAEmHFpcl6xQv
0WJV5iVVOq/ZgHKFpM8Ah0t6F3AT/unfSvASLVZFXlKlO5oKKBHxtyTLonyJZAn7/x0RnyoyY9ZZ
Hk5pVeLvc3c02yn/AWBNRHw4Ij4UEV8vOF/WYXntzLUmA9dWrBdkv6/uN+mOZmfKXwi8FdgGfBH4
z4h4uOC8TZiHDbem/ieEPaPYekXeDHjwkirjVfhPAKcXeRlwLvA7wKaIeF2rF+wkB5TxW7UqaX/e
vTv5K2/lyn19LWZl4+9re3XiJ4Ah+R2Uh0hmtB/R6sWsd7jJwHqJv6/l0GyT1/kkTV7PBf4TuCIi
1hSctwlzDWVi6pvAvEKxlU125WBwE1e7FNrkJeljwBcj4vvjyVy3OKC0h1dotTLy97I4hTZ5RcRH
gGdIemd6sedKOrbVi1lv8hBMKyN/L8vHa3nZmDyk2Mqk9t2rrRzsfpPyaLbJ6/vAiSQ/w3tiuu/u
iHhZwfmbEDd5tY+HFFsZeOXgzih6lJfX8upztSVapk1zU4N1j1cOLjev5WUtq28Cq/1Il5u/rAie
Ad87WvmBrdcDZwACbuiF5Vfc5FUc/0iXdYJnwHdHR2bK9xoHlOJ5hrIVyd+v7iikD0XSTkk7ch47
Je0Yf3atKtz8ZUXwSK7eVHgNRdIC4BKS4LU8Ii6qOz4V+DxwMvAT4NyIeFDSdJIl808BPhcR78+8
5hbg+cDPSQYKnBERP8m5tmsoHeDmL2snj+Tqvk6t5dUSSZOAS4EzgTnAYkkn1CU7D9gWEbNJAs/F
6f4ngY8CH2xw+sURcWJEnJQXTKxz/Bv11k4eydW7Cg0owDxgJCI2RMQuYAWwsC7NQuCydPtK4HSA
iPhZRNwG/KLBuYvOu7XIEyBtvDySqxqKLpRnABszzzel+3LTRMQeYHva3DWWz0q6Q9JH25JTm7C8
36j3zwrbWPJ+rte/Ad+bpnQ7Azmaabd7W0RsTSdYflnS2yMidymYwcHBvdsDAwMMDAy0JZOWL/sb
9atWHdgE5hE6Vi9vouz8+f6udNLQ0BBDQ0MTPk+hnfKS5gODEbEgfX4BENmOeUnXpWlulzQZ2BoR
R2SOLwFOznbK112j4XF3yndX7S/PNWuSpotrr036WbwEvsG+pednzUoGc9S+J66VdN94O+WLrqEM
A8dJmgVsBRYBi+vSXA0sAW4HzgFuzjnP3jeWBp3DI+JRSQcBZwGln2TZj2pNYB4BZvU8kquaOjVs
+BPsGzb8cUnLgOGIuEbSwcDlJItPPgosioj16WsfAKYBU4HtJDP1HwRWkgTDySTLwHwgryriGkp5
5E1QmzPHP9rVb2q1kp/+FN7wBk9YLCvPlM/hgFIeec1frrH0l2yt5IR08sB997mZq4wcUHI4oJRL
dgn8e+7Zv8Zy3XXw9Ke7tlI12Z/ozfs/P/RQN3OVkQNKDgeU8srWWI4/Ptm3bp1rK1WS10/izvfe
4ICSwwGl3Go1lieecHt6FTXqN/NKweVXyqVXzEZTm7Ny6qleYLJKRlvYMftDbVY9rqFYKXiByWrw
cOBqcJNXDgeU3uPhxb3Jw4GrxU1eVgl5v6/i9cDKLbsW15/+aTIk2As79ifXUKx0PLy4/DwcuNrc
5JXDAaX3eXhx+Xg4cPU5oORwQKmGRsOLXVvprNH6STwcuFocUHI4oFSLayvd42VT+osDSg4HlOoZ
bTKkR4O1l/tJ+pdHeVlfGG0yZHY02JYtnhw5HrVJiVu27P951k9SnDfPExTtQK6hWM9qNBpsyhQ4
5hhYv97NYa3INmvVPj/3k/QnN3nlcEDpH9n+lVmz8gtDN4fly+tsrwXlDRvcT9KPHFByOKD0l/rl
Wxr99kq//xRxtm8EGne2e9mU/uWAksMBpX812xzWL8Gl/vfba5/N3/7tgUOx3dluDig5HFAMGjeH
5fW1QDWaxhrVQur7Rq67Dj74QU9KtP05oORwQLGavOaw+r6WWuHaq01jzdRC8vpGwJ3ttr/SDhuW
tEDSOkn3S1qac3yqpBWSRiStkjQz3T9d0s2Sdkr6ZN1rTpJ0d3rOS4p+D9b7asONjzoqKURXroRv
fnP/obARScG6e3fy72mn7b8oZW1IbZmGIucN8z3ttH3vY80akPa9zzlzkve9cuW+2oh/o8TapdAa
iqRJwP3A6cAWYBhYFBHrMmneA7w0Is6XdC7w5ohYJOnpwCuAucDciHh/5jW3A++LiGFJ1wKfiIgb
cq7vGoqNKtvXAo1Hio1We4Fim8nqm6/GaspyLcQmarw1FCKisAcwH7gu8/wCYGldmuuBU9PtycAj
dceXAJ/MPD8SWJN5vgj4dIPrh1krduyIWLUqYvPmiJe/POKgg5J/v/71iClTIiD597jjkn/nzk0e
U6Yk6XbsSB633Tb6du1aYx2r5SPvWnl5quV38+bkfdTOZ9aKtOxsucyfMoEg1owZwMbM803AvEZp
ImKPpO2SpkfEtlHOuanunDPalF/rc7XmH0j+ss/WXubMObD2snZt0qRUa15avXpfTaY2DHfduv23
syv01qerP5atedRfq9aUlR0enR3me9RRHfzgzKDwgDIerVezRjE4OLh3e2BggIGBgXae3iosG1xg
X4DJduzXFqmszd3I9sNkA0B9MPja1/LT1R9bv35f81X9tebN2z/oOYjYeA0NDTE0NDTh8xTdhzIf
GIyIBenzC0iqUhdl0lyXprld0mRga0QckTm+BDg50j4USUcCt0TEi9Pni4DTIuI9OdePIt+f9a/6
vpe8fphsAKgPBtnfEBntWH3NI3st94VYUUo5bDgNEPeRdMpvBVYDiyNibSbN+SSd7uenweHsiFiU
Ob4E+NWI+KPMvu8A7yfp5P8aSR/L9TnXd0Cxjhst2GSDQaN09cccOKzTShlQIBk2DHyCZIjy8oj4
uKRlwHBEXCPpYOBy4ETgUZJRYOvT1z4ATAOmAtuBMyJinaSTgX8FDgGujYg/bnBtBxQzsxaVNqB0
kwOKmVnrSjux0czM+oMDipmZtYUDipmZtYUDipmZtYUDipmZtYUDipmZtYUDipmZtYUDipmZtYUD
ipmZtYUDipmZtYUDipmZtYUDipmZtYUDipmZtYUDipmZtYUDipmZtYUDipmZtYUDipmZtYUDipmZ
tUXhAUXSAknrJN0vaWnO8amSVkgakbRK0szMsY+k+9dKOiOzf72kuyTdKWl10e/BzMzGNqXIk0ua
BFwKnA5sAYYlfTUi1mWSnQdsi4jZks4FLgYWSXoJ8FbgxcDRwE2SZqc/Ev8UMBARjxWZfzMza17R
NZR5wEhEbIiIXcAKYGFdmoXAZen2lcBr0+03ASsiYndErAdG0vMBCDfXmZmVStGF8gxgY+b5pnRf
bpqI2AM8Lml6zms3Z14bwA2ShiW9q4iMm5lZawpt8honNZHm1yNiq6TnAl+XtDYivlV0xszMrLGi
A8pmYGZ7se6FAAAFhklEQVTm+dHpvqxNwAuALZImA4dFxDZJm9P9B7w2Iram/z4i6SskTWG5AWVw
cHDv9sDAAAMDAxN4O2Zm1TM0NMTQ0NCEz6Okj7sYaYC4j6RTfiuwGlgcEWszac4H5kbE+ZIWAWdH
RK1T/gvAqSRNXV8HZgNPAyZFxBOSDgVuBJZFxI05148i35+ZWRVJIiKaaS3aT6E1lIjYI+l9JIX+
JGB5RKyVtAwYjohrgOXA5ZJGgEeBRelr10i6AlgD7ALOj4iQ9DzgK5Iizf8X8oKJmZl1VqE1lG5z
DcXMrHXjraF46K2ZmbWFA4qZmbWFA4qZmbWFA4qZmbWFA4qZmbWFA4qZmbWFA4qZmbWFA4qZmbWF
A4qZmbWFA4qZmbWFA4qZmbWFA4qZmbWFA4qZmbWFA4qZmbWFA4qZmbWFA4qZmbWFA4qZmbWFA4qZ
mbWFA4qZmbVF4QFF0gJJ6yTdL2lpzvGpklZIGpG0StLMzLGPpPvXSjqj2XOamVnnFRpQJE0CLgXO
BOYAiyWdUJfsPGBbRMwGLgEuTl/7EuCtwIuBNwD/qEQz57Q2Gxoa6nYWKsWfZ3v58yyHomso84CR
iNgQEbuAFcDCujQLgcvS7SuB16bbbwJWRMTuiFgPjKTna+ac1ma+YdvLn2d7+fMsh6IDygxgY+b5
pnRfbpqI2AM8Lml6zms3p/uaOaeZmXVYGTvl1e0MmJlZ66YUfP7NwMzM86PTfVmbgBcAWyRNBg6L
iG2SNqf761+rJs65l+T41C7Lli3rdhYqxZ9ne/nz7L6iA8owcJykWcBWYBGwuC7N1cAS4HbgHODm
dP9VwBck/T1Jk9ZxwGqSWtVY5wQgIhxNzMw6pNCAEhF7JL0PuJEkECyPiLWSlgHDEXENsBy4XNII
8ChJgCAi1ki6AlgD7ALOj4gAcs9Z5PswM7OxKSmjzczMJqaMnfLjJuktku6RtEfSSaOk88TIMUh6
lqQbJd0n6QZJz2yQbo+kOyTdKem/Op3PspvIxF7bXxOf5RJJP06/j3dI+v1u5LNXSFou6WFJd4+S
5pPpd/P7kl4x1jkrFVCAHwBvBr7ZKIEnRjbtAuCmiDiepF/rIw3S/TQiToqIEyPi7M5lr/wmMrHX
9tfCfbsi/T6eFBGf7Wgme8/nSD7PXJLeALww/W7+T+CfxjphpQJKRNwXESOMPvTYEyObk51wehnQ
KFh44ENj45nYe3oH89dLmr1v/X1sUkR8C3hslCQLgc+naW8HninpeaOds1IBpUmeGNmcIyLiYYCI
eAg4okG6gyWtlnSbJAfm/Y1nYu/2dGKv7a/Z+/a30+aZKyQd3ZmsVVajyeUNFT1suO0kfR3IRkkB
AfxFRFzdnVz1plE+y4/mJG80emNWRGyVdCxws6S7I+KBNme1n/gv7PG7Cvj3iNgl6d0kNT/X+Dqo
5wJKRLx+gqdoZrJlXxjts0w7654XEQ9LOhL4cYNzbE3/fUDSEHAi4ICSGPfE3g7lr5eM+VlGRLb5
5l9wf9RENZpc3lCVm7wa/aW3d7KlpKkk816u6ly2esZVwO+l20uAr9YnkHR4+hki6TnAK0nmDVmi
me9abWIv7D+x1/Y35meZ/uFTsxB/F5shGpeVVwHvAJA0H9heawZvKCIq8yDpON4I/JxkFv116f7n
A9dk0i0A7iNZwfiCbue7jA9gOnBT+jndCBye7j8Z+H/p9q8BdwN3AncBv9ftfJftkfddA5YBZ6Xb
BwNXpMe/AxzT7TyX9dHEZ/nXwD3p9/EbwIu6necyP4B/B7YAvwAeBN5JMprr3Zk0lwL/nd7fJ411
Tk9sNDOztqhyk5eZmXWQA4qZmbWFA4qZmbWFA4qZmbWFA4qZmbWFA4qZmbWFA4qZmbWFA4qZmbXF
/wegN+xEoM28FwAAAABJRU5ErkJggg==
"
>
</div>
</div>
</div>
</div>
</div>
<div class="cell border-box-sizing text_cell rendered">
<div class="prompt input_prompt">
</div>
<div class="inner_cell">
<div class="text_cell_render border-box-sizing rendered_html">
<p>Note that when the correct weights are used, the points nearest zero have a great influence on the fitted model in their vicinity - their uncertainties are so small (and therefore their weights are so high) they're going to pull the best-fit line strongly towards them. With uniform weights, the endpoints have more pull just because of their influence on the slope.</p>
</div>
</div>
</div>
<div class="cell border-box-sizing text_cell rendered">
<div class="prompt input_prompt">
</div>
<div class="inner_cell">
<div class="text_cell_render border-box-sizing rendered_html">
<h3 id="robust-standard-errors">Robust standard errors</h3>
<p>Fortunately, linear least-squares is a very well-studied problem, and there's a nice matrix formula that is supposed to produce a very similar sort of calculation. Fundamentally what it tries to do is use the residual from each data point to estimate the uncertainty that data point should have had. Of course, this does not converge to a sensible estimate of the individual uncertainties, since no matter how many data points you have you have just as many uncertainties to estimate. But the uncertainties on your fit parameters are a finite set of modest size, and there's hope that an appropriate calculation could converge to those uncertainties. There are a few slightly different formulae for calculating this. I'm using the Mackinnon and White estimator, as described in <a href="http://www.indiana.edu/~jslsoc/files_research/testing_tests/hccm/99TAS.pdf">Long and Ervin 1999</a> (as HC3). This is <span class="math display">\[(A^TA)^{-1} A^T \text{diag}\left(\frac{r_i^2}{1-h_i^2}\right)_i (A^TA)^{-1},
\]</span> where <span class="math inline">\(r_i\)</span> is the <span class="math inline">\(i\)</span>th residual, and <span class="math inline">\(h_i\)</span> is the leverage of the <span class="math inline">\(i\)</span>th point.</p>
</div>
</div>
</div>
<div class="cell border-box-sizing code_cell rendered">
<div class="input">
<div class="prompt input_prompt">In [17]:</div>
<div class="inner_cell">
<div class="input_area">
<div class="highlight"><pre><span class="n">res</span> <span class="o">=</span> <span class="n">linear_least_squares</span><span class="p">(</span><span class="n">Acols</span><span class="p">,</span> <span class="n">ys</span><span class="p">,</span> <span class="n">claimed_uncerts</span><span class="p">)</span>
<span class="k">print</span> <span class="s2">"Robust standard errors:"</span>
<span class="k">print</span> <span class="n">res</span><span class="o">.</span><span class="n">cov_robust</span>
<span class="k">print</span> <span class="s1">'"True", approximated:'</span>
<span class="k">print</span> <span class="n">cov_het</span>
</pre></div>
</div>
</div>
</div>
<div class="output_wrapper">
<div class="output">
<div class="output_area"><div class="prompt"></div>
<div class="output_subarea output_stream output_stdout output_text">
<pre>
Robust standard errors:
[[ 3.78233363e-04 -7.43917638e-05]
[ -7.43917638e-05 1.79503265e-03]]
"True", approximated:
[[ 3.38980301e-03 -5.14196812e-05]
[ -5.14196812e-05 1.82065550e-02]]
</pre>
</div>
</div>
</div>
</div>
</div>
<div class="cell border-box-sizing code_cell rendered">
<div class="input">
<div class="prompt input_prompt">In [17]:</div>
<div class="inner_cell">
<div class="input_area">
<div class="highlight"><pre>
</pre></div>
</div>
</div>
</div>
</div>
</div>
</div>Unknownnoreply@blogger.com0tag:blogger.com,1999:blog-1369432396898204613.post-19537947832570623362016-07-11T15:54:00.000+02:002016-07-11T15:54:02.759+02:00EurovisionSince I now live in Europe, I find myself drawn in to the Eurovision song contest. On the one hand, it's an international contest even more fraught with under-the-surface (and on-the-surface) <a href="http://www.thewire.com/entertainment/2014/05/your-guide-to-the-politics-of-the-eurovision-song-contest/361707/">politics</a> than the <a href="http://www.theguardian.com/politics/politicspast/page/0,9067,892902,00.html">Olympics</a>; one is drawn to root for one's team. But on the other hand, it is a festival of massively-produced pop music. I mean, the live performances use the very highest technology available — <a href="http://www.eurovision.tv/page/news?id=stage_design_for_2016_revealed">massive LED matrices</a>, <a href="https://animationport.wordpress.com/2016/05/17/augmented-reality-projection-mapping-eurovision-song-contest-russia/">realtime projection mapping</a>, <a href="https://youtu.be/e94dst20C9Y">elaborate robotic sets</a>, even <a href="https://youtu.be/43h9aeyGL0E">goofy two-wheeled contraptions</a> — and the music is processed to within an inch of its life. Which is actually kind of fun to watch. Of course, all the songs are horrible earworms. So be warned; the video below is lovely and apropos for this blog, but it will be stuck in your head:<br>
<iframe allowfullscreen="" frameborder="0" height="315" src="https://www.youtube.com/embed/yBrADG8lWFY" width="560"></iframe><br>
Oh, and of course, it's astonishingly gay — the hosts opened with a string of jokes about all the gay fans, and I saw at least as many rainbow flags as I did flags from any nation (should I make that any <em>geographical</em> <a href="https://en.wikipedia.org/wiki/Queer_Nation">nation</a>?).<br>
<br>
Below the jump I'll post a few others that I particularly liked (that is, they're hopelessly stuck in my head).<br>
<br>
<a href="http://lighthouseinthesky.blogspot.com/2016/07/eurovision.html#more">Read more »</a>Unknownnoreply@blogger.com0tag:blogger.com,1999:blog-1369432396898204613.post-81694639946083388642016-07-10T18:20:00.000+02:002016-07-10T18:20:03.048+02:00In defense of Duck and CoverNow that the Cold War is over, people will occasionally look back at the horrifying situation — two opposing superpowers threatening the world with annihilation — view the threat as over, and laugh at some of the things people said and did back then. Now I know that this is partly because laughter is the only way to deal with such a nightmare, but some of the specifics they choose are just wrong. Take for example this John Oliver clip about nuclear weapons:
<iframe allowfullscreen="" frameborder="0" height="315" src="https://www.youtube.com/embed/1Y1ya-yF35g" width="560"></iframe><br>
He argues convincingly that there is still a threat, and that we need to do something about it. But he opens by mocking "Duck and Cover", which is actually <a href="https://en.wikipedia.org/wiki/Duck_and_cover#Efficacy_during_a_nuclear_explosion">misguided</a>.
<br>
<br>
<a href="http://lighthouseinthesky.blogspot.com/2016/07/in-defense-of-duck-and-cover.html#more">Read more »</a>Unknownnoreply@blogger.com0tag:blogger.com,1999:blog-1369432396898204613.post-55406503381864858212016-07-09T20:00:00.000+02:002016-07-09T22:43:25.478+02:00I for one welcome my new robot overlords, or at least housematesI recently gave in to temptation and obtained a 3D printer. Combined with my <a href="https://www.youtube.com/watch?v=ewdbilSWjaM">robot vacuum cleaner</a>, that means I am now outnumbered, and will be even more so when Coral (the 3D printer) has printed itself enough upgrades and starts printing parts for the upcoming Windeybot. Expect, therefore, many robot-oriented posts in the next little while. I'm still considering getting a cat to swing the balance a little towards organic life, but until then probably no cute cat posts. I'll see what I can do otherwise.<br>
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjYfW7MCW6l92CNTfxHCMMgU0hPeE69OIFWNmsFNMv09SwlM0G4xwkHO_g-a8fKfA9ACDqZs0Pmep9zME5w5CRqUHOCGLefLJl17JwmL4wZc9tEZdJcw4amFUltDMnLTqDXqT0tWhy7q_U/s1600/IMG_20160708_234446.jpg" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="640" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjYfW7MCW6l92CNTfxHCMMgU0hPeE69OIFWNmsFNMv09SwlM0G4xwkHO_g-a8fKfA9ACDqZs0Pmep9zME5w5CRqUHOCGLefLJl17JwmL4wZc9tEZdJcw4amFUltDMnLTqDXqT0tWhy7q_U/s640/IMG_20160708_234446.jpg" width="480"></a></div>
<br>
<br>
The 3D printer is from <a href="https://builda3dprinter.eu/">Builda3DPrinter</a>, a Dutch company (actually just one guy living in Hardinxveld) that sells <a href="http://builda3dprinte.wpengine.com/3d-printer-kits/kossel-mini/">kits</a> and <a href="https://octofiber.com/3d-printing-filament/pla.html">supplies</a> for <a href="http://www.parallemic.org/Reviews/Review002.html">delta</a> <a href="https://www.youtube.com/watch?v=PPyet-kDWMQ">printers</a>. There's lots more to say about the robot and how it works, but as a hobby, 3D printing is an interesting mix of getting the machine to work well and coming up with cool things to print. I'll just say a bit about it below.<br>
<br>
<br>
<a href="http://lighthouseinthesky.blogspot.com/2016/07/i-for-one-wolcome-my-new-robot.html#more">Read more »</a>Unknownnoreply@blogger.com0tag:blogger.com,1999:blog-1369432396898204613.post-84109353526285959202015-10-03T01:20:00.001+02:002015-10-03T01:20:04.684+02:00Numerical derivatives in parallelMy <a href="http://lighthouseinthesky.blogspot.nl/2015/10/numerical-coroutines.html">last post</a> was about a fairly esoteric new language feature in python: coroutines. I wasn't reading up on them because they're whiz-bang cool (okay, not <i>only</i> because they're whiz-bang cool) but because I have an actual problem that they might help with. For the <a href="http://lighthouseinthesky.blogspot.nl/2014/02/new-pet-psr-j033717.html">triple system</a>, I find I need to do very slow numerical optimization problems. That is, the objective function takes seconds to evaluate. It can't be parallelized internally, and numerical minimization is largely a sequential operation: evaluate the function, figure out which way is downhill, go there, repeat. So this code currently uses just one core, and is almost as slow as the <a href="http://dan.iel.fm/emcee/current/">highly parallel Markov-Chain Monte Carlo code</a> that does hundreds of times as much work. So obviously it would be nice to find some parallelism. I have a way to do that: parallelize gradient evaluations.<br>
<br>
The other reason for this post is to make a point about code robustness. Getting numerical derivatives right is a notoriously <a href="https://en.wikipedia.org/wiki/Numerical_differentiation#Practical_considerations_using_floating_point_arithmetic">difficult</a> process, and getting them both efficient and right is even harder. The package I have been using, <a href="http://seal.web.cern.ch/seal/snapshot/work-packages/mathlibs/minuit/">MINUIT</a>, has been used by the particle physicists for decades, and they have hammered out all sorts of corner cases that a naive implementation might miss. So I'm not going to use a naive implementation, not least because my problem is numerically difficult. This means the code is not simple, and I don't want to try to restructure it manually to be more concurrent; I want the computer to handle the concurrency for me.<br>
<br>
<a href="http://lighthouseinthesky.blogspot.com/2015/10/numerical-derivatives-in-parallel.html#more">Read more »</a>Unknownnoreply@blogger.com2tag:blogger.com,1999:blog-1369432396898204613.post-5492065129700849782015-10-03T00:50:00.001+02:002015-10-03T01:25:12.012+02:00Numerical coroutinesPython 3.5 has <a href="http://makina-corpus.com/blog/metier/2015/python-http-server-with-the-new-async-await-syntax">new syntax</a> for coroutines. Like many of python3's new features (<a href="http://www.diveintopython3.net/strings.html">unicode strings</a>? <a href="http://www.astronomie.nl/#!/home/">who</a> uses <a href="https://craq-astro.ca/">non-English</a> <a href="http://adsabs.harvard.edu/abs/1964AZh....41..282K">languages</a> in astronomy?), it is not immediately obvious that coroutines are of any use to people like me who mostly write relatively simple imperative scripts to wrangle <a href="http://lighthouseinthesky.blogspot.nl/2009/05/finding-pulsars-part-2.html">large amounts of data</a> or <a href="http://lighthouseinthesky.blogspot.nl/2014/10/parallel-ipython-notebooks.html">long computations</a>. But cores aren't getting any faster, so we're going to have to start writing parallel code if we're going to take advantage of new machines. Now, if you want to take a <a href="http://link.springer.com/chapter/10.1007/978-1-4612-1516-5_14">gigantic Fourier transform</a>, or do <a href="http://www.osti.gov/scitech/biblio/71638">algebra on enormous matrices</a>, you are going to need a witch design you a parallel algorithm. But our code is often full of easy opportunities for parallelism that we don't take advantage of because it's awkward. So I keep my eye out for ways to take advantage of this easy parallelism. I have found one, in the combination of coroutines and futures, and I'd like to describe it with a worked example. Maybe it'll sell you on some of these cool new language features!<br>
<br>
<a href="http://lighthouseinthesky.blogspot.com/2015/10/numerical-coroutines.html#more">Read more »</a>Unknownnoreply@blogger.com0tag:blogger.com,1999:blog-1369432396898204613.post-30066599279879349592014-10-01T21:20:00.001+02:002014-10-05T17:18:26.862+02:00Parallel ipython notebooks and the H test
<div class="text_cell_render border-box-sizing rendered_html">
<p>The CPU on the nice shiny new server I log in to is really not much faster than that of the ratty old laptop I have in front of me. The server has more memory, and more disk space, but CPU-wise it just has <em>more</em> not faster. For that matter even my laptop has two cores. So if I have some heavy-duty computing task, I'd better find a way to make it use multiple cores in parallel. Some tasks are just plain hard to parallelize (<a href="http://lighthouseinthesky.blogspot.co.uk/2014/02/new-pet-psr-j033717.html">solving ordinary differential equations</a>, for example), but it turns up fairly often that I'm doing something <a href="http://lighthouseinthesky.blogspot.co.uk/2009/05/finding-pulsars-part-2.html">embarrassingly parallel</a>: there's a for loop somewhere that just does the same thing to lots of different pieces of input. If only it were easy to hand each piece to a different core! Well, there are <a href="http://pymotw.com/2/multiprocessing/index.html">various</a> <a href="http://dan.iel.fm/emcee/current/user/advanced/#using-mpi-to-distribute-the-computations">tools</a> for doing this sort of thing, but most of them apply to scripts or programs that run non-interactively. It turns out that ipython offers <a href="http://ipython.org/ipython-doc/2/parallel/parallel_intro.html">tools for <em>interactive</em> parallel computing</a>. I'm going to explain how I use them, by working through a test problem, checking some statistics on a periodicity test (the H test).</p>
</div><a href="http://lighthouseinthesky.blogspot.com/2014/10/parallel-ipython-notebooks.html#more">Read more »</a>Unknownnoreply@blogger.com0tag:blogger.com,1999:blog-1369432396898204613.post-20512107146256306622014-09-29T11:50:00.000+02:002014-09-29T11:50:43.455+02:00Remote ipython notebooks<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjS5cgqhWAFLDIljaQgUuVY6ZtLb3sIRUnv0UyrwK6zA4tEfDvIxsfOE5mQ_Ab3GAXMrjhlZ-4lHpi_zztVrxiJJWPxO3P8oNV846rjYnBYsQoYDFyOTwF7FAqNqr3_KnLfn4fxCFPUuuI/s1600/Screenshot+from+2014-09-29+09:52:36.png" imageanchor="1" style="clear: right; float: right; margin-bottom: 1em; margin-left: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjS5cgqhWAFLDIljaQgUuVY6ZtLb3sIRUnv0UyrwK6zA4tEfDvIxsfOE5mQ_Ab3GAXMrjhlZ-4lHpi_zztVrxiJJWPxO3P8oNV846rjYnBYsQoYDFyOTwF7FAqNqr3_KnLfn4fxCFPUuuI/s1600/Screenshot+from+2014-09-29+09:52:36.png" height="111" width="200"></a></div>
I use ipython notebooks a lot. For interactive tinkering, for calculations with units, for interactive control of long-running processes, for easy parallelization (details to come), and to generate plots for all my papers. When I can, I run the notebooks on a server rather than my rather wimpy laptop; apart from the obvious advantages of more cores, more disk space, and more RAM, this also means I can leave things running for ages. It's all very convenient, but it requires some setup. As with the ssh configuration, I've found myself explaining the setup to several different people recently, so I thought I'd put it here. (Still working on that cute-kittens post.)<br>
<br>
<a href="http://lighthouseinthesky.blogspot.com/2014/09/remote-ipython-notebooks.html#more">Read more »</a>Unknownnoreply@blogger.com0tag:blogger.com,1999:blog-1369432396898204613.post-57254568553850346892014-09-25T20:15:00.001+02:002014-09-29T11:34:11.931+02:00SSH incantationsI use ssh all the time, for shell access on remote machines, to transfer files, to carry VNC sessions, and to use different institutional affiliations to get at paywalled things. Over the years I have developed a particular way of setting it up that minimizes the pain involved in doing all of the above things. Several times in the past months I have tried to explain the whole setup to someone, so I think that's a sign it's time to post it here. I'll try to post something about, I don't know, fluffy kittens, boating on the Thames, <a href="http://pipeline.corante.com/archives/2014/09/23/seaborgium_hexacarbonyl.php">seaborgium hexacarbonyl</a>, or something soon for those of you who for some reason expected an interesting blog.<br>
<br>
<a href="http://lighthouseinthesky.blogspot.com/2014/09/ssh-incantations.html#more">Read more »</a>Unknownnoreply@blogger.com2tag:blogger.com,1999:blog-1369432396898204613.post-3788407962062626402014-03-03T21:00:00.000+01:002014-03-04T23:06:55.786+01:00The Fortress of the North<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEh_drrn0BHV_o8nBzSmO_iC1jel9W6AbGiK0Du-Pvud-kTuMwRZxz6zAGaRsSIpd9I09zB8ddcTtu_siY3vJ_joJd7vZaQeZ1zoKrpQ1yfGNXxrHjbbW2ht3-qfzW5H6joJghQwKC2aD3I/s1600/IMG_3991.JPG" imageanchor="1" style="clear: right; float: right; margin-bottom: 1em; margin-left: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEh_drrn0BHV_o8nBzSmO_iC1jel9W6AbGiK0Du-Pvud-kTuMwRZxz6zAGaRsSIpd9I09zB8ddcTtu_siY3vJ_joJd7vZaQeZ1zoKrpQ1yfGNXxrHjbbW2ht3-qfzW5H6joJghQwKC2aD3I/s1600/IMG_3991.JPG" height="240" width="320"></a>I live in Groningen. It's a city of about 200,000, about 50,000 of whom are students at <a href="http://www.hanze.nl/home/International">the</a> <a href="http://www.rug.nl/?lang=en">universities</a>. So I sort of thought it might resemble <a href="http://en.wikipedia.org/wiki/Kitchener,_Ontario">Kitchener</a>-<a href="http://en.wikipedia.org/wiki/Waterloo,_Ontario">Waterloo</a>, the college town where I did my undergraduate degree, but in fact I've never been anywhere quite like this before.<br>
<br>
<a href="http://lighthouseinthesky.blogspot.com/2014/03/the-fortress-of-north.html#more">Read more »</a>Unknownnoreply@blogger.com1tag:blogger.com,1999:blog-1369432396898204613.post-37979312772578316882014-03-02T18:00:00.000+01:002014-03-02T18:00:02.339+01:00The problem with "beginner's mind"I should say up front that I know nothing about authentic Zen Buddhism. But one of the concepts that has made it into <a href="http://www.amazon.com/Zen-Art-Motorcycle-Maintenance-Inquiry/dp/0060589469">popular culture</a> that I rather like is "<a href="http://en.wikipedia.org/wiki/Shoshin">beginner's mind</a>": an <a href="http://www.goodreads.com/quotes/306037-none-of-our-men-are-experts-we-have-most-unfortunately">expert</a> knows what she thinks about things, while a beginner sees everything anew and judges it on what it is. It sounds kind of charming, and like a recipe for creativity. But applying it too literally can get you in trouble.<br>
<br>
Concretely, <a href="http://www.astron.nl/">ASTRON</a> (where I <a href="http://lighthouseinthesky.blogspot.nl/2013/08/new-job-new-country.html">work</a>) arranged for a <a href="http://www.atcomputing.nl/">company</a> to come in and teach programming courses to anyone who wanted them. They taught three: Introductory Python, Numeric Python, and Advanced Python. Now, I think that's a brilliant idea: we astronomers all spend most of our time writing code at one level or another, and nobody seems to have bothered to teach us how to do it well. So that ASTRON (actually <a href="http://www.nwo.nl/en">NWO</a> I think) cares is a really good sign. So of course I wanted to participate. I have lots to learn about writing good programs. I don't think I'm an expert programmer, or even a Python expert, but that latter is partly because I try not to think of myself as an expert in anything. I knew the Introductory Python course would not be productive — I <a href="http://lighthouseinthesky.blogspot.nl/2011/09/review-ipython-notebooks.html">write</a> python code every day. And the Numerical Python, again, was a good idea, but having written, for example, the <a href="http://docs.scipy.org/doc/numpy/reference/generated/numpy.reshape.html">reshape</a> <a href="https://github.com/numpy/numpy/blob/688b243552cc98d8e1be791eeda3fba81bfbe807/numpy/core/src/multiarray/shape.c">function</a> in numpy, and the <a href="http://docs.scipy.org/doc/scipy/reference/spatial.html">spatial</a> <a href="https://github.com/scipy/scipy/blob/master/scipy/spatial/ckdtree.pyx">module</a> in scipy, I figured that was probably not going to be too productive either. But the Advanced Python course sounded promising. And I wanted to show my support for the whole idea. So I signed up.<br>
<br>
Well, the course was the last three days, and it was a pretty good course, but not at all what I needed. Which should not have been much of a surprise: when I stopped to do the math, I realize that I wrote my first python program eighteen years ago. (<a href="http://en.wikipedia.org/wiki/History_of_Python#Version_release_dates">Good God.</a>) Still, it was interesting to see how they ran the course, and to think about how I would run one (because if I end up somewhere like McGill that has basically nothing for physics grad students, I will run one, official or not).<br>
<br>
<a href="http://lighthouseinthesky.blogspot.com/2014/03/the-problem-with-beginners-mind.html#more">Read more »</a>Unknownnoreply@blogger.com1tag:blogger.com,1999:blog-1369432396898204613.post-47317452888198082842014-03-01T18:00:00.000+01:002014-03-02T14:11:46.934+01:00Visualizing the new petI <a href="http://lighthouseinthesky.blogspot.nl/2014/02/new-pet-psr-j033717.html">recently</a> wrote about a new object I am studying: a <a href="http://lighthouseinthesky.blogspot.nl/2009/05/where-do-millisecond-pulsars-come-from.html">millisecond pulsar</a> with two white dwarf companions. There is lots more I want to say about it, but I think it would be nice to show what it looks like, or at least, to show a video I made trying to make visible what's going on:<br>
<div class="separator" style="clear: both; text-align: center;">
<object class="BLOGGER-youtube-video" classid="clsid:D27CDB6E-AE6D-11cf-96B8-444553540000" codebase="http://download.macromedia.com/pub/shockwave/cabs/flash/swflash.cab#version=6,0,40,0" data-thumbnail-src="https://i1.ytimg.com/vi/oDgfqq_W_uM/0.jpg" height="266" width="320"><param name="movie" value="https://www.youtube.com/v/oDgfqq_W_uM?version=3&f=user_uploads&c=google-webdrive-0&app=youtube_gdata"><param name="bgcolor" value="#FFFFFF"><param name="allowFullScreen" value="true"><embed width="320" height="266" src="https://www.youtube.com/v/oDgfqq_W_uM?version=3&f=user_uploads&c=google-webdrive-0&app=youtube_gdata" type="application/x-shockwave-flash" allowfullscreen="true"></object></div>
<i>Edited to note that Blogger's YouTube embedding is distinctly flaky; video is <a href="https://www.youtube.com/watch?v=oDgfqq_W_uM">here</a>. </i><br>
<i><br></i>
This video shows the orbital motions in the triple system. The orbits are drawn to scale, showing the actual motions of the two stars (red and yellow) and the pulsar (white). The first ten seconds are played relatively slowly, showing the motion around the inner orbit, then we speed up to see the motion around the outer orbit. For a sense of the time scale, an "MJD" is a modified Julian day, so a single day long. The larger left panel shows all three bodies, with trails marking the motion of the outer companion and the center of mass of the inner system. The inset in the top right zooms in on the inner system, showing the pulsar and the companion, with trails marking their orbits. The dots that appear on the orbits mark moments when we have observations of the system, color-coded by telescope; it should be clear that we have quite thorough coverage of both orbits. Each measurement tells us the distance to the pulsar to within a kilometer, so that we can measure the tiny deviations of these orbits from perfect Keplerian ellipses, allowing us to reconstruct the orbit.<br>
<br>
There's a little more to it than that.<br>
<br>
<a href="http://lighthouseinthesky.blogspot.com/2014/03/visualizing-new-pet.html#more">Read more »</a>Unknownnoreply@blogger.com0tag:blogger.com,1999:blog-1369432396898204613.post-60715786029522933422014-02-18T21:00:00.000+01:002014-02-18T21:00:00.412+01:00New pet: PSR J0337+17<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjuRHXzwiOYhyphenhyphend1bAntHW9pvniTg-OF-h3X1H_wK_SV-mT4ry93usupEbd4N9eFazfwqucZnHF9aPvJU1xxIz-GmiyHeQvcn4nidJPcL0x8bPUmf0yZtCrbMybofTNfrtOV4596M0QGQI4/s1600/0337+17-profile.png" imageanchor="1" style="clear: right; float: right; margin-bottom: 1em; margin-left: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjuRHXzwiOYhyphenhyphend1bAntHW9pvniTg-OF-h3X1H_wK_SV-mT4ry93usupEbd4N9eFazfwqucZnHF9aPvJU1xxIz-GmiyHeQvcn4nidJPcL0x8bPUmf0yZtCrbMybofTNfrtOV4596M0QGQI4/s1600/0337+17-profile.png" height="256" width="320"></a>I did my PhD thesis on <a href="http://lighthouseinthesky.blogspot.nl/2009/05/missing-link.html">PSR J1023+0038</a>, a millisecond pulsar that is at a fascinating point in its evolution. (In fact there have been developments since the thesis was submitted; more about that later.) But during a moment of procrastination, I got involved with a new and fascinating system. The name, unedifying as usual, is PSR J0337+17, and it is unique in that the pulsar has not just one white dwarf companion but two.<br>
<br>
<br>
<br>
<a href="http://lighthouseinthesky.blogspot.com/2014/02/new-pet-psr-j033717.html#more">Read more »</a>Unknownnoreply@blogger.com1tag:blogger.com,1999:blog-1369432396898204613.post-29551807088815954772014-02-17T21:00:00.000+01:002014-02-17T21:00:00.112+01:00Flywheel energy storageIn the quest for something better to run our cars on than gasoline, one of the <a href="http://www.youtube.com/watch?v=VowB4MGQqQQ">proposals</a> is <a href="http://www.youtube.com/watch?v=eCtlfj4kMJs">flywheels</a>. In fact, for a while there were flywheel-powered buses <a href="http://en.wikipedia.org/wiki/Gyrobus#Early_commercial_service">running in Switzerland and Belgium</a>. On one level, it makes a lot of sense: you're storing energy as mechanical motion, and we're pretty good at <a href="http://www.youtube.com/watch?v=7kUPoQu-jjE">transmitting</a> mechanical motion from place to place. On another level it scares the living daylights out of me: a car in motion uses tens of kilowatts, so the car must be able to store hundreds of kilowatt-hours. If you let all those loose at once bad things will happen: 100 g of TNT going off releases about a hundred kilowatt-hours. Fortunately it's <a href="http://jalopnik.com/why-cars-explode-into-fireballs-and-why-they-usually-do-560552028">hard</a> to get gasoline to do this, but a flywheel is just <a href="http://www.youtube.com/watch?v=y7gKbk0jyyM">itching</a> to <a href="http://www.chem.purdue.edu/chemsafety/newsandstories/centrifugedamages.htm">dump</a> all its energy. Batteries are a little <a href="http://www.youtube.com/watch?v=k9mcNvOGKtI">scary</a> too, to be honest. But anyway, that's all a digression: I want to talk about some really staggering examples of flywheel energy storage: pulsars and black holes.<br>
<br>
<a href="http://lighthouseinthesky.blogspot.com/2014/02/flywheel-energy-storage.html#more">Read more »</a>Unknownnoreply@blogger.com0tag:blogger.com,1999:blog-1369432396898204613.post-78727375732808745202013-08-28T21:00:00.000+02:002013-08-28T21:00:05.553+02:00New job, new countryI just started a new job, and in the way of academic life, I had to move. I put some things in storage, hopped on a plane with two suitcases, and hey presto, I live in the Netherlands now. The new job is with <a href="http://www.astron.nl/">ASTRON</a>, the Dutch institute for radio astronomy, and it's going really well. Living in a new country actually takes more getting used to.<br>
<br>
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgCXMa3khK3Lw_jsn3dElnJupJNWnjROFGZ6gW4r28iak5-W7BhS-5uMMKnSJn_4oX37yLc2LWnLHR0QKIwyjveND7IZnPab93OsRqaus0LQoz1m_eLN846YeoAZnMU5jdQ1yai7xCnIdI/s1600/IMG_3548.JPG" imageanchor="1" style="clear: right; float: right; margin-bottom: 1em; margin-left: 1em;"><img border="0" height="240" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgCXMa3khK3Lw_jsn3dElnJupJNWnjROFGZ6gW4r28iak5-W7BhS-5uMMKnSJn_4oX37yLc2LWnLHR0QKIwyjveND7IZnPab93OsRqaus0LQoz1m_eLN846YeoAZnMU5jdQ1yai7xCnIdI/s320/IMG_3548.JPG" width="320"></a>For one thing, in spite of the fact that I live in arguably the <a href="http://en.wikipedia.org/wiki/Pedestrian_zone#Pedestrianised_centres">most pedestrianized</a> <a href="http://en.wikipedia.org/wiki/Groningen">city</a> in Europe, and in a densely-populated country with good transit systems and cities built when rapid transit meant corn-fed horses, I find myself tempted to carpool to work. See, the thing is, ASTRON was built in the middle of a <a href="http://en.wikipedia.org/wiki/Dwingelderveld_National_Park">national park</a> in order to minimize interference with the <a href="http://en.wikipedia.org/wiki/Dwingeloo_Radio_Observatory">radio telescope</a> that was also being built. Since the fifties, though, radio astronomy has made a few strides, and the telescope is now used chiefly in education. The <a href="http://www.lofar.org/">telescopes</a> <a href="http://www.astron.nl/radio-observatory/astronomers/wsrt-astronomers">we</a> <a href="http://images.nrao.edu/Telescopes/GBT">do</a> <a href="http://www.naic.edu/public/the_telescope.htm">use</a> are all <a href="http://fermi.gsfc.nasa.gov/">off-site</a>, so we're in the middle of the park for historical reasons. It means we're twelve kilometers from the <a href="http://en.wikipedia.org/wiki/Beilen_railway_station">nearest train station</a>, and even the bus drops you off in the village of <a href="http://en.wikipedia.org/wiki/Lhee">Lhee</a>, a kilometer or so away from work. So a ride to work is tempting. Fortunately, this is the Netherlands, so I've been biking.<br>
<br>
<a href="http://lighthouseinthesky.blogspot.com/2013/08/new-job-new-country.html#more">Read more »</a>Unknownnoreply@blogger.com3tag:blogger.com,1999:blog-1369432396898204613.post-14110518946832056282013-02-26T02:00:00.000+01:002013-02-26T02:00:01.845+01:00Differentiating the solutions to differential equations<div class="text_cell_render border-box-sizing rendered_html">
A kind of problem that turns up quite often in physics and elsewhere is to find the solution of an ordinary differential equation, given some initial conditions. That is, you have a system with some state $x$ represented as a vector of real numbers, you have an initial state $x_0$, and you have a rule describing the evolution of the state:
$$
\frac{dx}{dt} = F(x,t)
$$
And your goal is to find $x(t)$. This standard problem has some standard solution techniques, some quite advanced - <a href="http://en.wikipedia.org/wiki/Dormand%E2%80%93Prince_method">Runge-Kutta methods</a>, <a href="http://en.wikipedia.org/wiki/Symplectic_integrator">symplectic methods</a>, <a href="http://adsabs.harvard.edu/abs/1992PASJ...44..141M">Hermite integrators</a>. A few are implemented <a href="http://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.ode.html#scipy.integrate.ode">in scipy</a>. But it sometimes happens that solving this problem is only part of a <a href="http://en.wikipedia.org/wiki/Shooting_method">more complicated process</a>, say of <a href="http://lighthouseinthesky.blogspot.ca/2010/12/how-tempo2-does-its-fitting.html">fitting</a>, where it would be nice to have the derivatives of the solution with respect to the various initial conditions. It turns out this isn't too hard to work out, usually.<br>
<br>
</div><a href="http://lighthouseinthesky.blogspot.com/2013/02/differentiating-solutions-to.html#more">Read more »</a>Unknownnoreply@blogger.com0tag:blogger.com,1999:blog-1369432396898204613.post-54526717019642621952013-02-02T00:03:00.000+01:002013-02-07T23:12:31.142+01:00Weighted Poisson Uncertainties<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhJBVxoQL557mG65L_ZzP0hJYBcNeqBP1G1yeHlTXDbDX8ggLzjU-UxvAxyoPKpxAUVuxg8VpRanW9zaZ86tYFWMdUUiHO-AHpI8orvC3-CNlOZsvXmEUOlnjQrrlhsj9x6Q_dSfwDQr9s/s1600/grpulse.png" imageanchor="1" style="clear: right; float: right; margin-bottom: 1em; margin-left: 1em;"><img border="0" height="218" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhJBVxoQL557mG65L_ZzP0hJYBcNeqBP1G1yeHlTXDbDX8ggLzjU-UxvAxyoPKpxAUVuxg8VpRanW9zaZ86tYFWMdUUiHO-AHpI8orvC3-CNlOZsvXmEUOlnjQrrlhsj9x6Q_dSfwDQr9s/s320/grpulse.png" width="320"></a>I recently ran across a rather awkward mathematical problem. I'm trying to make a histogram of photon arrival phases, complete with an uncertainty on the number of photons in each bin. Normally this is done by just taking the square root of the number of photons, which is at least approximately right based on Poisson statistics. But in my problem — data from the Fermi space telescope, which is not very good at localizing low-energy gamma rays — the photons are weighted: for each photon I have a probability that it really came from the source. So the values in the histogram should be the total probability. But what should the uncertainty be? The short version is: the square root of the sum of the squares of the weights.<br>
<br>
<b>ETA:</b> This is in the literature, without justification as far as I can tell. See below.<br>
<a href="http://lighthouseinthesky.blogspot.com/2013/02/wieghted-poisson-uncertainties.html#more">Read more »</a>Unknownnoreply@blogger.com0tag:blogger.com,1999:blog-1369432396898204613.post-83341454630677255642012-12-21T06:01:00.000+01:002012-12-21T06:01:21.021+01:00git tear-hairI use <a href="http://git-scm.com/">git</a>. It's fast, it's convenient, it gives you <a href="https://github.com/">github</a>, and mostly it works. But several times now I've managed to bollix a git repository, been faced with impenetrable messages, and been unable to continue using git unless I sort them out. This latest time I had enough stubbornness to figure out how to un-bollix my repository, and I wanted to record it for posterity (specifically, for people whose problem-solving technique involves Googling for the specific error message). The error I received was:<br>
<pre>error: object file .git/objects/86/5d2dffe9a3d72917934ed9693c7167efb6d8d5 is empty
fatal: loose object 865d2dffe9a3d72917934ed9693c7167efb6d8d5 (stored in .git/objects/86/5d2dffe9a3d72917934ed9693c7167efb6d8d5) is corrupt
</pre>
Read on for how it happened, my understanding of what it means, and how I fixed it. Or skip this even-more-technical-than-usual post; I'll try to post something with <a href="http://www.crosscreekutila.com/photos/Acropora-cervicornis.jpg">corals</a> or <a href="http://travel-to-honduras.com/photo/1507-85/Three_Roatan_Kittens_in_a_Row.htm">kittens</a> soon.<br>
<a href="http://lighthouseinthesky.blogspot.com/2012/12/git-tear-hair.html#more">Read more »</a>Unknownnoreply@blogger.com0tag:blogger.com,1999:blog-1369432396898204613.post-32269477269967494162012-12-05T05:08:00.001+01:002012-12-05T05:10:11.893+01:00Justice<table cellpadding="0" cellspacing="0" class="tr-caption-container" style="float: left; margin-right: 1em; text-align: left;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjmCtzC_EXDmvpuJXxVjOY4Zsoi26LuNPYFaPS9VIXQ8h14TjlOr-_M3KipPZziueQ_7Njq1ULEwGZkajMhWl7TyIrb0L4wS-5PDnGjY2MmNW3K7XvDlG48YT9Bven6Mb8kfuSLBgej2iE/s1600/justice.jpg" imageanchor="1" style="clear: left; margin-bottom: 1em; margin-left: auto; margin-right: auto;"><img border="0" height="320" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjmCtzC_EXDmvpuJXxVjOY4Zsoi26LuNPYFaPS9VIXQ8h14TjlOr-_M3KipPZziueQ_7Njq1ULEwGZkajMhWl7TyIrb0L4wS-5PDnGjY2MmNW3K7XvDlG48YT9Bven6Mb8kfuSLBgej2iE/s320/justice.jpg" width="320"></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">Available <a href="http://www.amazon.com/Justice-Whats-Right-Thing-Do/dp/0374532508">from Amazon</a></td></tr>
</tbody></table>
I recently read (listened to? What's the right way to say that for an audio book?) "<a href="http://books.google.ca/books/about/Justice.html?id=BrdNDG7TTUEC">Justice: What's the right thing to do?</a>" by <a href="http://www.gov.harvard.edu/people/faculty/michael-sandel">Michael J. Sandel</a>. It's an interesting book, though rather heavier going than Nate Silver's <a href="http://lighthouseinthesky.blogspot.ca/2012/11/prediction.html">breezy book on prediction</a>. The basic question it addresses is, how do you decide what's right? Not just what's legal — you need some basis for deciding which laws are just. He talks through several approaches, including some ideas of his own. While he talks about several currently-debated topics, he refrains from stating his position on any of them, instead pointing out what the conflicting ideas of justice are that motivate the two sides.<br>
<br>
<a href="http://lighthouseinthesky.blogspot.com/2012/12/justice.html#more">Read more »</a>Unknownnoreply@blogger.com0tag:blogger.com,1999:blog-1369432396898204613.post-84169314556352738052012-11-16T02:26:00.003+01:002012-11-16T02:26:39.743+01:00Talking about radio astronomy<div class="separator" style="clear: both; text-align: center;">
<br /></div>
<div class="separator" style="clear: both; text-align: center;">
<iframe allowfullscreen='allowfullscreen' webkitallowfullscreen='webkitallowfullscreen' mozallowfullscreen='mozallowfullscreen' width='320' height='266' src='https://www.youtube.com/embed/XvwNCNd85IA?feature=player_embedded' frameborder='0'></iframe></div>
McGill has a modest <a href="http://www.physics.mcgill.ca/~aarchiba/observatory.html">observatory</a> on the roof, with a 14-inch (optical) telescope. It was installed quite some time ago, then left to molder for a few years. This annoyed me, just on principle, so I talked my way into rehabilitating it as best I could, with my radio astronomer's skills. Even with it working again, it was just used by the occasional grad student who wanted to take her friends up to see Jupiter, or <a href="http://lighthouseinthesky.blogspot.ca/2011/11/lucky-imaging.html">the surface of the moon</a>, or whatever. Fortunately, we recently acquired an enthusiastic and organized post-doc who pushed hard and set up a <a href="http://www.astro.physics.mcgill.ca/outreach.php">public outreach program</a>. Typical <a href="http://www.astro.physics.mcgill.ca/astronight.php">nights</a> start with a public lecture, then we take folks outside to two portable telescopes and the bigger one on the roof to look at the stars. Of course, we can't predict the weather, so we have a few demos we can do inside - comet-making, liquid nitrogen, a muon detector, but really the appeal is going out and looking up at the sky. But especially in the summer, we need to entertain the public until it gets dark. Since it's once a month, we have a perennial need for speakers. So I volunteered, to give a one-hour talk on radio astronomy for the general public. I had never given a one-hour talk before, and radio astronomy doesn't produce as many pretty pictures as I would like, but I think it came out pretty well. And we had a videographer who recorded the whole thing, so if you're curious, <a href="http://www.youtube.com/watch?v=XvwNCNd85IA">here</a> it is.Unknownnoreply@blogger.com0tag:blogger.com,1999:blog-1369432396898204613.post-29827077698569534372012-11-12T09:45:00.000+01:002012-11-12T11:25:14.737+01:00Prediction<table cellpadding="0" cellspacing="0" class="tr-caption-container" style="float: right; margin-left: 1em; text-align: right;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEh3cuG2p5SPY34jmSqHTwLVPzGl4J1VAuLTvIqhQwIgun1uubA9nCkLkpFCDCN1hKkl6By5mWZP-R-PCdkAt3Wx6ZISq_PNA7UVXY45xdkaOgd0YCMMvT6mcNcCIU6ogcOEu25FyMFu4bE/s1600/51KLaBetObL._SS500_.jpg" imageanchor="1" style="clear: right; margin-bottom: 1em; margin-left: auto; margin-right: auto;"><img border="0" height="320" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEh3cuG2p5SPY34jmSqHTwLVPzGl4J1VAuLTvIqhQwIgun1uubA9nCkLkpFCDCN1hKkl6By5mWZP-R-PCdkAt3Wx6ZISq_PNA7UVXY45xdkaOgd0YCMMvT6mcNcCIU6ogcOEu25FyMFu4bE/s320/51KLaBetObL._SS500_.jpg" width="320"></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">The book can be bought <a href="http://www.amazon.com/dp/159420411X">from Amazon</a></td></tr>
</tbody></table>
I recently finished reading <a href="http://en.wikipedia.org/wiki/Nate_Silver">Nate Silver</a>'s book "The Signal and The Noise: Why Most Predictions Fail But Some Don't". Nate Silver is currently somewhat famous because his election predictions - he runs the blog <a href="http://fivethirtyeight.blogs.nytimes.com/">fivethirtyeight</a> for the New York Times - his election predictions were basically spot on. The book, timed to come out shortly before the election (in what I think was a sort of gamble to capitalize on a successful prediction) is about prediction, generally defined. Structured as a sort of collection of case studies, the book talks about all the ways prediction can go wrong, and about what to do to try to get it right. There's no nice one-line answer, but I think he gives some pretty good advice.<br>
<br>
<a href="http://lighthouseinthesky.blogspot.com/2012/11/prediction.html#more">Read more »</a>Unknownnoreply@blogger.com5tag:blogger.com,1999:blog-1369432396898204613.post-67679163092035964002012-09-29T18:00:00.000+02:002012-09-29T18:00:05.653+02:00The Maya-Muon experimentThis week's Friday colloquium talk really had me feeling like I was in a science fiction show. The <a href="http://www.ph.utexas.edu/~schwitte/Schwitters.htm">speaker</a> is <a href="http://www.hep.utexas.edu/mayamuon/aboutus/">building machines</a> to use cosmic-ray muons to see the interiors of still-buried Mayan ruins. Does this not sound like a line of technobabble from Stargate SG-1? But it's nearly feasible.<div>
</div><a href="http://lighthouseinthesky.blogspot.com/2012/09/the-maya-muon-experiment.html#more">Read more »</a>Unknownnoreply@blogger.com1tag:blogger.com,1999:blog-1369432396898204613.post-76641069290413067372012-07-04T09:08:00.000+02:002012-07-04T11:00:18.772+02:00Higgs liveblog<table cellpadding="0" cellspacing="0" class="tr-caption-container" style="float: right; text-align: right;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEicFfmhcVYz1bsDA1Ol9XdAkcpo2t0Kn9cy7-C49MTstuBZauRgyCM8wwyLv9nzgmPm5QWmFHyrHAL4jBYF0wynks0Jg2NlIJfvwAQmmUq7bJK6JXxMiJ2Pkj6tXKibZ99MrsUMpJCTx30/s1600/cand-3.png" imageanchor="1" style="clear: right; margin-bottom: 1em; margin-left: auto; margin-right: auto;"><img border="0" height="145" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEicFfmhcVYz1bsDA1Ol9XdAkcpo2t0Kn9cy7-C49MTstuBZauRgyCM8wwyLv9nzgmPm5QWmFHyrHAL4jBYF0wynks0Jg2NlIJfvwAQmmUq7bJK6JXxMiJ2Pkj6tXKibZ99MrsUMpJCTx30/s200/cand-3.png" width="200"></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">A candidate Higgs event; see below.</td></tr>
</tbody></table>
The Higgs detection (?) announcement is being <a href="http://webcast.web.cern.ch/webcast/play_higgs.html">broadcast live</a>. So I thought I'd blog my comments as it unfolds.<br>
<br>
<a href="http://lighthouseinthesky.blogspot.com/2012/07/higgs-liveblog.html#more">Read more »</a>Unknownnoreply@blogger.com0tag:blogger.com,1999:blog-1369432396898204613.post-42020379041095494112012-07-03T02:00:00.000+02:002012-07-03T02:00:03.724+02:00FireworksIn honor of Canada Day (and for my neighbours to the south, the fourth of July) here's a video:<br>
<div class="separator" style="clear: both; text-align: center;">
<object class="BLOGGER-youtube-video" classid="clsid:D27CDB6E-AE6D-11cf-96B8-444553540000" codebase="http://download.macromedia.com/pub/shockwave/cabs/flash/swflash.cab#version=6,0,40,0" data-thumbnail-src="http://0.gvt0.com/vi/u-Xp_-_gLTA/0.jpg" height="266" width="320"><param name="movie" value="http://www.youtube.com/v/u-Xp_-_gLTA&fs=1&source=uds">
<param name="bgcolor" value="#FFFFFF">
<param name="allowFullScreen" value="true">
<embed width="320" height="266" src="http://www.youtube.com/v/u-Xp_-_gLTA&fs=1&source=uds" type="application/x-shockwave-flash" allowfullscreen="true"></object></div>
While this might look like a meteor, and in fact it is asteroidal material falling to Earth, it's actually some extremely expensive fireworks. It's the <a href="http://en.wikipedia.org/wiki/Hayabusa_(spacecraft)">Hayabusa spacecraft</a> returning to Earth after visiting asteroid <a href="http://en.wikipedia.org/wiki/25143_Itokawa">25143 Itokawa</a>, and the material it brought back was <a href="http://hayabusa.seti.org/">safely encapsulated</a>.<br>
<br>
<a href="http://lighthouseinthesky.blogspot.com/2012/07/fireworks.html#more">Read more »</a>Unknownnoreply@blogger.com0tag:blogger.com,1999:blog-1369432396898204613.post-84371432826728951572012-06-30T02:00:00.000+02:002012-06-30T02:00:02.470+02:00The Too-wide WebFor some reason, computer monitors have been getting wider and wider for years. This puzzles me, since like most people, when I'm working, I tend to use tall narrow documents, both to read and to edit. Sometimes I can arrange things so that I have two panels on the screen, which restores them to a more sensible shape, but the Web is a problem. Web pages seem to have begun adapting to wide monitors by adding wider and wider margins, often filled with ads and/or navigation materials. For me, these margins are often too wide for me to use a two-panel setup (it's just a laptop) but often they leave so much width when used full-screen that the text is tiresomely long. Typesetters have a <a href="http://www.maxdesign.com.au/articles/em/">rule of thumb</a> that you shouldn't put more than about twelve words on a line because it's hard to read. Fortunately, I found a Chrome hack that lets me solve the problem.<br>
<table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjQgL86LVN0Hpvsz-d3J_RnssnetDQMrvmz1HvtOBH2KOLNgL5WpCF98cd1rdAonldT9MD4Mn52Jse5pv2wMU9AqaPasQADKmyEIhjfofG77HloKGp7xMzyJajnt6yHSP__YG2UL7JwZAA/s1600/Screenshot+from+2012-06-29+12:29:53.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="179" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjQgL86LVN0Hpvsz-d3J_RnssnetDQMrvmz1HvtOBH2KOLNgL5WpCF98cd1rdAonldT9MD4Mn52Jse5pv2wMU9AqaPasQADKmyEIhjfofG77HloKGp7xMzyJajnt6yHSP__YG2UL7JwZAA/s320/Screenshot+from+2012-06-29+12:29:53.png" width="320"></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">The <a href="http://reprap.org/wiki/Planning_and_preparing_for_a_build">RepRap wiki</a> is too wide.</td></tr>
</tbody></table>
<br>
<br>
<a href="http://lighthouseinthesky.blogspot.com/2012/06/too-wide-web.html#more">Read more »</a>Unknownnoreply@blogger.com4