Microsoft Word – Problem Set #3 Solutions.doc
MA 568 – Statistical Analysis of Point Process Data
Solutions to Problem Set 3
% Question 1
load Hipp_data.mat
s1Ind = find(spikes);
s2Ind = find(spikes2);
plot(xN,yN,xN(s1Ind),yN(s1Ind),’.’,
xN(s2Ind),yN(s2Ind),’.’);
axis tight square;
These two neurons are both tuned to the spatial
location of the animal, and both have similar
tuning properties in that their place fields nearly
completely overlap. The center of both place
fields occurs when both the x and y-coordinates
are slightly negative.
% Question 2
[b,dev1,stats] = glmfit([xN yN],spikes,’poisson’);
stats.p
[x_new,y_new]=meshgrid(-1:.1:1);
lambda = exp(b(1) + b(2)*x_new + b(3)*y_new);
h_mesh = mesh(x_new,y_new,lambda,’AlphaData’,0);
All of the p values are highly significant (p<10-26). The exponentiated linear model captures the fact that the firing rate increases as we move from the origin to negative values for x and y. However, the firing rate continues to increase in this direction, even though the number of spikes falls off as x and y become more negative. In order to capture this phenomenon, we will need to include quadratic terms of x and y. % Question 3 X = [xN yN xN.^2 yN.^2 xN.*yN]; [b,dev2,stats] = glmfit(X,spikes,'poisson'); stats.p >> stats.p = [0 0 0 0 0 0.6196]
All of the parameters except the one corresponding to the xN.*yN term are significant.
This suggests that a more parsimonious model could be constructed by removing this
term from the model.
-0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8
-1
-0.8
-0.6
-0.4
-0.2
0
0.2
0.4
0.6
0.8
1
page 2: MA 568 – Solutions to Problem Set 3
lambda = exp(b(1) + b(2)*x_new + b(3)*y_new + b(4)*x_new.^2 +
b(5)*y_new.^2 + b(6)*x_new.*y_new );
h_mesh = mesh(x_new,y_new,lambda,’AlphaData’,0);
AIC1 = dev1 + 2*3
AIC2 = dev2 + 2*6
>> AIC1 = 11394
>> AIC2 = 9153.2
This Gaussian place field model is able to more
accurately describe the firing properties of each
of these neurons. Each place field has a center
where firing rate is maximal, and falls off
evenly when moving in any direction away. Even with the extraneous xN.*yN term, the
Gaussian model significantly lowers the AIC value.
% Question 4
for p = 0:20,
[b,dev,stats] = glmfit([X spikes_hist(:,1:p)],spikes,’poisson’);
AIC(p+1) = dev+2*length(b);
end;
plot(0:20,AIC);
The AIC is minimized by a history component that
goes back 8 ms. The optimal model from this
model class would include the constant term, the 4
quadratic terms with significant parameters, and 8
history terms.
% Question 5
[b,dev,stats] = glmfit([X spikes_hist(:,1:8) spikes2_hist],
spikes, ‘poisson’);
stats.p
None of the parameter estimates related to the history of the second neuron are
significant. This suggests that there is no functional interaction from neuron 2 to neuron
1, once position related effects have been modeled. The optimal model remains the one
with a constant term, the 4 quadratic terms with significant parameters, 8 history terms,
and no interaction terms.
% Question 6
[b,dev,stats] = glmfit([spikes_hist(:,1:8) spikes2_hist], spikes,
‘poisson’);
stats.p
When we eliminate the position related terms from the model, all of the interaction
parameters become significant. This is not surprising, since any spiking from neuron 2 is
likely to correspond to periods when the animal is near the center of place field 2, and
0 2 4 6 8 10 12 14 16 18 20
8000
8200
8400
8600
8800
9000
9200
A
IC
V
al
ue
Number of History Terms
page 3: MA 568 – Solutions to Problem Set 3
therefore, near the center of the place field of neuron 1. This model uses information
about position that is implicit in the firing of neuron 2 to predict the spiking of neuron 1.
When we add the position term, this information is already present and no significant
interaction terms are observed.
% Question 7
[b,dev,stats] = glmfit([X spikes_hist(:,1:8)],spikes,’poisson’);
lambda_ML = exp([ones(size(T)) X spikes_hist(:,1:8)]*b);
lambdaInt = cumsum(lambda_ML);
Z = diff([0; lambdaInt(s1Ind)]);
U = expcdf(Z,1);
n = length(U);
M = [1:n]/n;
plot(sort(U),M,M,M+1.36/sqrt(n),’k:’,M,M-1.36/sqrt(n),’k:’);
axis([0 1 0 1]);
The most parsimonious model does not completely
pass the KS test. This suggests that there still may be
other statistical structure in the spiking activity that
cannot be captured within this model class.
However, the KS statistic is relatively small
suggesting that this model is able to explain most of
the observed structure in the data.
0 0.2 0.4 0.6 0.8 1
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Empirical CDF
M
od
el
C
D
F