x=[1:10]'
p=[0.5,1]
y=polyval(p,x)
e(1:5,1)=0.1;
e(6:10,1)=1;
for i=1:100
yd=y+e.*randn(size(y));
[pf(i,:),pfe(i,:)]=matmin('chisq',p,[],'polyval',yd,e,x);
end
subplot(2,1,1); hfill(pf(:,1));
subplot(2,1,2); hfill(pf(:,2));
[std(pf(:,1)) pfe(1,1)]
[std(pf(:,2)) pfe(1,2)]
Note that pfe(:,1) is a column of equal numbers - the parameter errors depend only on the data errors and hence are equal for each sim realization. Also note that the fir values scatter around the true values and that the spread (std) of the scatter matches the parameter error estimate.
The example we will use is counts in a histogram. If we have n counts in a bin people often take the sqrt(n) as the error on the bin and do the fit chisq style as x,n,sqrt(n) as x,y,e. But what is the chance to get a count of 1? 1+/-1? No if we see even a single event clearly the probability to get zero is zero!
What is the probability to get 1 event if the model says the expectation value in that bin is 0.6 events?
poisspdf(1,0.5)
Matlab has an excellent set of probability functions. Look at the Statistics toolbox help under "By Category".
[x,n]=hfill(randn(1,1000),100,-3,3);
hplot(x,n);
[p,pe]=matmin('chisq',[30,0,1],[],'gauss',n,sqrt(n),x);
hold on
plot(x,gauss(p,x),'r');
hold off
function L = logl(pars,func,n,varargin)
% L = logl(pars,func,n,x...)
%
% Calculate -2 times the log of the joint probability of a
% set of observations n at values x if they are Poisson deviates from
% mean values described by function func with parameters
% pars.
%
% eg: [x,n]=hfill(randn(1,1000),100,-3,3);
% L=logl([30,0,1],'gauss',n,x)
% [p,pe]=matmin('logl',[30,0,1],[],'gauss',n,x)
% get model value for each n
mu=feval(func,pars,varargin{:});
% use poisspdf to get probability of each n given model
% combine these to get L
return
Now we fit using logl:
[p,pe]=matmin('logl',[30,0,1],[],'gauss',n,x)
hold on
plot(x,gauss(p,x),'g');
hold off
In this case the results are similar but it can maka a big difference when there are many bins with zero contents - chisq has to to be made to ignore them or divide by zero occurs.
[ra dec]=textread('wmap_ptsrc_catalog_p2_3yr_v2.txt','%s %s %*[^\n]','commentstyle','shell')
textread reads the entire file in one shot. It is very powerful but rapidly becomes convoluted. Check the Matlab help and Write a function to return a structure s with elements of glon, glat and a nx5 array of fluxes s=wmap_read('wmap_ptsrc_catalog_p2_3yr_v2.txt').
s=wmap_read('wmap_ptsrc_catalog_p2_3yr_v2.txt')
plot(s.l,s.b,'.'); xlim([0,360]); ylim([-90,90]);
clf
axesm('MapProjection','Hammer','Grid','on');
plotm(s.b,s.l,'.')
Now let's plot the histogram of source fluxes.
hfill(s.s(:,1),100);
The form of this becomes clearer if we plot in log bins with a log y axis;
s.s(s.s==0)=NaN; [x,nl]=hfill(log10(s.s(:,1)),30); hplot(x,nl); set(gca,'YScale','log'); ylim([0.3,100]);
Now let's overplot the highest frequency:
[x,nh]=hfill(log10(s.s(:,end)),30,x(1),x(end)); hold on; hplot(x,nh,'r'); hold off; ylim([0.3,100]);
At low flux we get no hits because the experiment is not sensitive enough. If we try and fit down there without telling the fitter about the selection efficiency we will get a silly answer. We need to limit to the region above flux where efficiency appears to become 100%.
[x,nl]=hfill(log10(s.s(:,1)),30,0); hplot(x,nl); set(gca,'YScale','log'); ylim([0.3,100]);