%
% Solution to exercise 1, parts 1 and 2, homework 2, part 2
% CAAM 415, 10/23/08
%
rho = 40; %in units of spk/sec

mu = 1000/40; %mean ISI in msec

figure(1);

for i = 1:10
    %generate twice as many spikes as expected to be on the safe
    %side
    r = exprnd(mu,1,2*rho); 
                            
    r2 = cumsum(r); %cumulative sum of intervals gives spike times

    spks = round(r2(find(r2<1000))); %find the spike occuring in the first 1000 msec

    %plot a series of lines for it
    x_spks = [1; 1]*spks;
    y_spks = [zeros(1,length(spks));ones(1,length(spks))];
    h = subplot(10,1,i);
    line(x_spks,y_spks,'Color','b');
 
    set(h,'YLim',[0 1],'YTick', [], 'YTickLabel',[]);
    if (i<10)
        set(h,'XTick',[]);
        set(h,'XTickLabel',[]);
    end;
        
end;

xlabel('time (ms)');

%Theoretical mean and variance
for j = 1:20
    t_int = 50*j;
    n_theor(j) = t_int/mu;
    v_theor(j) = t_int/(2*mu) + (1/8)*(1-exp(-t_int/mu));
end;

figure(2);
plot(n_theor,v_theor,'r');
hold on;

%compute the spike times from the ISI distribution. 
%the mean of the gamma is mu_g*n_g, so work out the right value of mu_g
n_g = 2;
mu_g = mu/n_g;

clear n2;
for i = 1:1000
    %generate twice as many ISIs as expected to be on the safe side
    r = gamrnd(n_g,mu_g,1,2*rho);
    %compute spike times from the ISIs by summing all ISIs up to the
    %current spike
    r2 = cumsum(r); 
    
    %round to 1 msec resolution and subselect only spike times in the first
    %second
    spks = round(r2(find(r2<1000)));
    
    for j = 1:20
        t_int2 = 50*j;
        n2(j).vals(i) = length(find(spks<t_int2));
    end;
end;

for j = 1:20
    n2_mean(j) = mean(n2(j).vals);
    v2_mean(j) = var(n2(j).vals);
    t_int = 50*j;
end;

plot(n2_mean,v2_mean,'xb');

xlabel('mean number of spikes');
ylabel('variance');

%plot FF as a function of time window
t_v = (1:20)*50;
ff = v_theor./n_theor;

figure(3);
plot(t_v,ff);