%
% Solution to exercise 2, homework 1, part 2
% CAAM 415, 10/16/08
%
load spont_data;
x1 = spont_data(:,1);
dx = x1(2)-x1(1); %recover the sampling interval

y1 = spont_data(:,2);

y1p = y1./sum(y1); %normalize to get probabilities
mu1 = sum(x1.*y1p); %mean depolarization

%reconstruct a vector of observations appropriate for
%normfit
x1_obs = [];
for i = 1:length(x1)
    x1_obs = [ x1_obs x1(i)*ones(1,y1(i))]; 
end;

%estimate mean and standard deviation    
[muhat,sigmahat,muci,sigmaci] = normfit(x1_obs);

%sample the Gaussian at higher resolution
dx_f = 0.01;
x_f = 0.0:dx_f:0.8;
yn = normpdf(x_f,muhat,sigmahat);

%normalize so that the area under both curves is identical
%(area under yn with sampling step dx_f is 1 by construction, 
%i.e., sum(yn*dx_f) = 1)
yn_rs = yn*sum(y1*dx);

%plot everything
figure(1);
bar(x1,y1);
hold on;
plot(x_f,yn_rs,'r');
xlabel('spontaneous e.p.p.s. (mV)');
ylabel('number of observations');

%same as above
load stim_data

x2 = stim_data(:,1);
y2 = stim_data(:,2);
n_x2 = length(x2);

y2p = y2./sum(y2);
mu2 = sum(x2.*y2p); %mean depolarization

l_e = mu2/mu1; %quantal size, standard method
l_f = -log(y2(1)/sum(y2)); %method of failures

%relative error
l_err = abs(l_e - l_f)/(0.5*(l_e + l_f));

info_str = sprintf('relative error: %.3f',l_err);
disp(info_str);

figure(2)
bar(x2,y2);

%sample again at high resolution
x3 = 0:dx_f:3.2;
n_x3 = length(x3);

y3 = zeros(1,n_x3);
y3(1) = poisspdf(0,l_e); %discrete probability of failures

%cycle over the peaks and implement the formula in the lecture notes
for i=1:7
    pois_val = poisspdf(i,l_e);
    mu = i*muhat; %mean of the ith peak
    sigma = sqrt(i)*sigmahat; %standard deviation
    yi = normpdf(x3,mu,sigma);
    y3= y3+pois_val*yi; %add everything toghether
end;

%normalize so that the area under everything but the first element 
%is equal to 1 - the probability of the 1st element (failures)
%this is the correct probability density
y3(2:n_x3)= y3(2:n_x3)*(1-y3(1))/sum(y3(2:n_x3)*dx_f);

%correctly normalized probability density
y3n = zeros(1,n_x3);

%the fraction of failures is the total sum of events times the probability
%of failures
y3n(1) = sum(y2)*y3(1);

%normalize so that the area under the continuous part is equal to 
%the area under the sample histogram
y3n(2:n_x3) = y3(2:n_x3)*(sum(y2(2:n_x2)*dx))/(sum(y3(2:n_x3)*dx_f));

hold on;
plot(x3,y3n,'r');
xlabel('evoked e.p.p.s. (mV)');
ylabel('number of observations');