<center><h2>
Class 3
</center></h2>

<center><img SRC="../../../colline.gif" ></center>

<ul>

<li>Last week's homework was to write a function fit_calsrc to auto
fit all the channels:

<pre>
load ../../050517_calsrc1_edit.mat
plot(t,v)
csfv=fit_calsrc(t,v)
hist(csfv(:,4),100)
</pre>

<p>The way sinusoid is written the phase is measured in cycles.
We are interested in phase as angle so we wrap it into the
range 0 to 1 and mult by 360:

<pre>
hist(mod(csfv(:,4),1)*360,100);
xlim([0,360]);
</pre>

<p>We see four groups roughly spaced by 90 deg.
The bolometers grid angles in the focal plane are separated by 45 deg steps -
why do we see 90 deg seps in the output waveforms? (Hint: pol is a "headless
vector".)

<li>Let's find and fix the glitch - which channel is it in?

<pre>
plot(v)
[mv,j]=min(min(v))
plot(v(:,j))
[mv,i]=min(v(:,j))
v(i-4:i+15,j)=NaN;
plot(v(:,j))
</pre>

<p>lsqnonlin does not treat NaN's as missing values (although regress does),
so we must filter them out ahead of time inside the loop of fit_calsrc.

<li>What if we know the uncertainty on each data point (and they are different)
and we want to include that knowledge when we make the fit?

<p>Make some "data" with differing errors:

<pre>
x=[1:10]'
p=[0.5,1]
y=polyval(p,x)
plot(x,y,'.')
e(1:5,1)=0.1;
e(6:10,1)=1;
y=y+e.*randn(size(y));
errorbar(x,y,e,'.');
</pre>

<p>Fit it first ignoring our knowledge of the uncertainties and again including
this knowledge -
remember lsqnonlin squares and sums the vector of numbers we feed it so
we are just minimizing chi-squared in the usual manner.

<pre>
pf1=lsqnonlin(@(v) y-polyval(v,x),p)
pf2=lsqnonlin(@(v) (y-polyval(v,x))./e,p)
hold on
plot(x,polyval(pf1,x),'r')
plot(x,polyval(pf2,x),'g')
hold off
</pre>

<p>Run it a few times and we see the benefit of down weighting the
higher uncertainty points.

<li>What if we want to assign uncertainties to the parameters?
lscov is a generalized regress and actually a better way to fit a linear
model poly above - and it also easily gives us errors on the parameters.
lsqnonlin can output the Jacobian which allows to do the same thing but
I prefer to use the more general MINUIT package which has a manual
I can understand and allows one to solve maximum likelihood
as well as least squared problems.
Here is the <a href="minuit.pdf">MINUIT manual</a> - have a look
at it - it contains a lot of good practical stuff.
<p>

<li>MINUIT is fortran code which can be interfaced to Matlab using
some special tricks to form a "Matlab plugin" binary or "mex"
file.
First we will look at a simple c code plugin which "fill" numbers
into a hitorgram where you get to control the bin range (unlike
Matlab hist function).
As we start to get more complicated we need to make a Matlab
toolbox directory and put it on our path:

<pre>
cd
mkdir matlab
cd matlab
emacs startup.m
</pre>

<p>Add the lines "dbstop if error" and "addpath('home/pryke/matlab']);"

<li>Now we get the hfillc.c and makefile:

<pre>
cp ~pryke/matlab/hfillc.c .
cp ~pryke/matlab/makefile .
emacs hfillc.c makefile &
make
</pre> 

<p>Now in Matlab window type "rehash" and "which hfillc" - we see it on
path.
At the moment to have it run we need to exit matlab, type 
"setenv LD_PRELOAD /lib/libgcc_s.so.1" and restart.
(This is a silly complication and best not to worry about why it's needed right now.)

<p>It's always best to put an m file wrapper around a mex file to do
things like argument checking etc:

<pre>
cp ~pryke/matlab/hfill.m .
cp ~pryke/matlab/hplot.m .
</pre>

<p>"rehash" and "help hfill" tells us about how to use it.

<li>Try to write your own version [bc,n]=hfill_alt(x,nbin,low,high)
using native Matlab code (i.e. without resorting to C-code extension).
Don't loop over x, but I think you do need to loop over nbin.
Time it by doing:

<pre>
x=randn(1e6,1);
tic; [bc1,n1]=hfill(x,100,-3,3); toc
tic; [bc1,n1]=hfill_alt(x,100,-3,3); toc
</pre>

<p>Can you get it anywhere near as fast as the C-code verison?
Can you avoid a loop over bins? (I don't know how.)

<li>Now let's get the MINUIT interface to Matlab.
First we need to copy the source code and compile.
Let's put it in our home directory in matlab sub-dir:

<pre>
cp -rp ~pryke/matlab/minuit .
make
</pre>

<p>We also need the wrapper and the "figure of merit functions":

<pre>
cp ~pryke/matlab/matmin.m .
cp ~pryke/matlab/chisq.m .
cp ~pryke/matlab/logl.m .
cp ~pryke/gauss.m .
</pre>

<p>Now type "rehash" and "help matmin" to see the instructions.
Finally we are all set to use MINUIT.
Let's go back to the example above and check it gives same result
as lsqnonlin:

<pre>
x=[1:10]'
p=[0.5,1]
y=polyval(p,x)
e(1:5,1)=0.1; e(6:10,1)=1;
y=y+e.*randn(size(y));
errorbar(x,y,e,'.');
pf=lsqnonlin(@(v) (y-polyval(v,x))./e,p)
pf=matmin('chisq',p,[],'polyval',y,e,x)
</pre>

<p>Except now we can get the errors on the parameters:

<pre>
[pf,pfe]=matmin('chisq',p,[],'polyval',y,e,x)
</pre>

<p>Note that MINUIT prints out the correlation matrix of the parameters.

<li>Are these errors correct? We can find out by doing our first Monte-Carlo...
This week's homework is to write code which will generate 100 sets of dithered
y's, fit each set using matmin, record the pf and pfe outputs, histogram the
100 pf(:,1) and pf(:,2) values, and compare the std of each against the
predicted pfe values from MINUIT.
You can also look at the correlation matrix of the fit parameters if you like.
If you want to be super clever modify matminc.c to fetch and return each
time its estimate of the parameter correlation matrix (which is
printed on the screen but not returned as an output).

</ul>

<hr>