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),’.’,
1
0.8
0.6
0.4
xN(s2Ind),yN(s2Ind),’.’); 0.2
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
0
-0.2
-0.4
-0.6
-0.8
fields occurs when both the x and y-coordinates -1
are slightly negative.
-0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8
% 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.
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);
9200
9000
8800
The AIC is minimized by a history component that 8600
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. 8000
% 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
8400
8200
0 2 4 6 8 10 12 14 16 18 20 Number of History Terms
AIC Value
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]);
1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
0 0.2 0.4 0.6 0.8 1
Empirical CDF
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.
Model CDF