MA 568 – Statistical Analysis of Point Process Data Solutions to Problem Set #1
% Question 1
load Retinal_ISIs.mat;
hist(ISIs,50);
spiketimes = cumsum(ISIs);
N = length(spiketimes);
for i = 1:N,
line([spiketimes(i)
spiketimes(i)],[0 1]);
end;
1 0.5 0
200 180 160 140 120 100
80 60 40 20
0
0 20 40 60 80 100 120 140 160 180 200
ISI (ms)
0 0.5 1 1.5 2 2.5 3
4 x10
The histogram and spike train time series suggest that this neuron fires most of its spikes with an ISI between 5-40 ms. Features to note include an initial refractory period under 5 ms., a large number of spikes with ISIs around 4-8 ms., a large tail that includes ISIs out to 190 ms., and a small second mode around 90-100 ms. This suggests that a simple Poisson model will not be able to fully capture the structure observed in this data.
% Question 2
Qs = sort(ISIs);
[Qs(1) Qs(end/4) Qs(end/2) Qs(3*end/4) Qs(end)]
%ans=2 5 10 44190 boxplot(ISIs);
1
0 20 40 60 80 100 120 140 160 180
The 5-number summary and boxplot confirm our earlier observations about the distribution of ISIs. Most spikes occur at small ISIs between 5-40 ms., but there are a large number of long ISIs over 100 ms.
page 2: MA 568 – Solutions to Problem Set #1
% Question 3
spikes1 = hist(spiketimes,0:30000);
plot(spikes1);
hist(spikes1,2);
4 x 10
2.5
2
13 0.5
0
1 0.5 0 1 0.5 0
1.5
1 0.5 0
0
1000
1.6
2000
1.7
3000
1.8
4000
1.9
5000
2
6000
2.1
7000
2.2
0.8
2.3
0.9
2.4
1
2.5
1.1
2.6
1.2
2.7
1.3
2.8
1.4
2.9
1.5
4 1.5 x 10
4 x 10
3
4 x 10
1
0.5
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
spikes10 = hist(spiketimes,0:10:30000);
plot(spikes10);
hist(spikes10,4);
4 3.5 3 2.5 2
2500
2000
1500
1.5 1000 1
500
0.5
00 0 0.5 1 1.5 2 2.5 3
4 x 10
01234
spikes100 = hist(spiketimes,0:100:30000);
plot(spikes100);
hist(spikes100,10);
10 9 8 7 6 5 4 3 2 1 0
60
50
40
30
20
10
0
lam = N/30
SE = sqrt(N)/30;
CI = [lam-1.96*SE lam+1.96*SE]
0 1 2 3 4 5 6 7 8 9 10
0 0.5 1 1.5 2 2.5 3
4 x 10
% Question 4
lams = [0:1e-1:60];
L = N*log(lams*1e-3)-lams*30;
plot(lams,L);
-4000 -4500 -5000 -5500 -6000 -6500 -7000 -7500 -8000 -8500 -9000
% lam = 32.4000
% CI = 30.3631
34.4369
0 10 20 30 40 50 60
page 3: MA 568 – Solutions to Problem Set #1
1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
% Question 5
w = [0:200];
for i = 1:length(w),
Femp(i) = sum(Qs<=w(i))/N;
end;
plot(w,Femp,w,1-exp(-w/lam))
Empirical CDF Model CDF
0 20 40 60 80 100 120 140 160 180 200
1 0.9 0.8 0.7 0.6 0.5 plot(Femp,1-exp(-w/lam)); 0.4 0.3 0.2 0.1 0
KSStat = max(abs(1-exp(-w/lam)-Femp))
% KSStat = 0.2314
axis([0 1 0 1]);
200
180
160
140
120 Qmodel = expinv(1/N:1/N:1,lam); 100
80 60 40 20
0
0 50 100 150 200
% Question 7
FF1 = var(spikes1)/mean(spikes1)
CI1 = gaminv([.025 .975], length(spikes1)/2, 2/length(spikes1))
% FF1 = 0.9676
% CI1 = 0.9841 1.0161
FF10 = var(spikes10)/mean(spikes10)
CI10 = gaminv([.025 .975], length(spikes10)/2, 2/length(spikes10))
% FF10 = 1.1271
% CI10 = 0.9500 1.0512
FF100 = var(spikes100)/mean(spikes100)
CI100 = gaminv([.025 .975], length(spikes100)/2, 2/length(spikes100))
% FF100 = 1.3411
% CI100 = 0.8464 1.1662
% Question 6
plot(Qs,Qmodel)
axis([0 200 0 200])
0 0.2 0.4 0.6 0.8 1
page 4: MA 568 – Solutions to Problem Set #1
% Question 8
Z = ISIs-mean(ISIs);
plot(-(N-1):(N-1),xcorr(Z,Z,'coef'),'.')
line([-N N],[1.96/sqrt(N) 1.96/sqrt(N)]);
line([-N N],[-1.96/sqrt(N) -1.96/sqrt(N)]);
1
0.8
0.6
0.4
0.2
0
-200 -150 -100 -50
0 50
100 150
200 250
At small lags, the autocorrelation function falls outside the 95% confidence bounds more often than expected by chance alone for small lags, suggesting that nearby spikes have ISIs that are not independent.
9. From the above analysis, we can conclude that under a simple Poisson model the firing rate of this neuron is likely somewhere between 30-34 Hz. However, this Poisson model fails to account for many aspects of the spiking activity including its refractoriness, bursting, and long ISIs. Our analysis of the Fano Factor suggests that there is less variability in the 1 ms. increments and more in the 10 ms and 100 ms increments than would be expected from a Poisson process. The autocorrelation analysis suggests a dependence structure on past spiking activity. Therefore, describing the structure of this data will require history dependent point process models.