Gardner Lab

Psychology Department

Neurosciences Institute

Stanford University


日本語   

Complete code for channel encoding model tutorial

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Make simulated data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
% some variables
iNeuron = 1;
orientations = 0:179;
k = 10;
% loop over each neuron tuning function
for orientPreference = 0:2:179  
  % compute the neural response as a Von Mises function
  %Note the 2 here which makes it so that our 0 - 180 orientation
  % space gets mapped to all 360 degrees
  neuralResponse(iNeuron,:) = exp(k*cos(2*pi*(orientations-orientPreference)/180));
  % normalize to a height of 1
  neuralResponse(iNeuron,:) = neuralResponse(iNeuron,:) / max(neuralResponse(iNeuron,:));
  % update counter
  iNeuron = iNeuron + 1;
end
 
% plot the response of neuron 45
figure;
plot(orientations,neuralResponse(45,:));
xlabel('Orientation');
ylabel('Channel response (normalized units to 1)');
 
% make a random weighting of neurons on to each voxel
nNeurons = size(neuralResponse,1);
nVoxels = 50;
neuronToVoxelWeights = rand(nNeurons,nVoxels);
 
% make stimulus array
nStimuli = 8;
% evenly space stimuli
stimuli = 0:180/(nStimuli):179;
% number of repeats
nRepeats = 20;
stimuli = repmat(stimuli,1,nRepeats);
 
% round and make a column array
stimuli = round(stimuli(:))+1;
 
% compute the voxelResponse
nTrials = nStimuli * nRepeats;
for iTrial = 1:nTrials
  % get the neural response to this stimulus, by indexing the correct column of the neuralResponse matrix
  thisNeuralResponse = neuralResponse(:,stimuli(iTrial));
  % multiply this by the neuronToVoxelWeights to get the voxel response on this trial. Note that you need
  % to get the matrix dimensions right, so transpose is needed on thisNeuralResponse
  voxelResponse(iTrial,:) = thisNeuralResponse' * neuronToVoxelWeights;
end
 
% plot the voxelResponse for the 7th trial
figure;
plot(voxelResponse(7,:));
xlabel('Voxel (number)');
ylabel('Voxel response (fake measurement units)');
 
% plot another trial voxel response
figure;
plot(voxelResponse(7,:),'b-.');
hold on
plot(voxelResponse(7+nStimuli,:),'r-o');
xlabel('Voxel (number)');
ylabel('Voxel response (fake measurement units)');
 
% add noise to the voxel responses
noiseStandardDeviation = 0.05;
% normalize response 
voxelResponse = voxelResponse / mean(voxelResponse(:));
% add gaussian noise
voxelResponse = voxelResponse + noiseStandardDeviation*randn(size(voxelResponse));
 
% check the voxelResponses
figure;
stim1 = 7;
stim2 = 3;
subplot(1,3,1);
plot(voxelResponse(stim1,:),'b-.');
hold on
plot(voxelResponse(stim1+nStimuli,:),'r-o');
xlabel('Voxel (number)');
ylabel('Voxel response (fake measurement units)');
 
subplot(1,3,2);
plot(voxelResponse(stim1,:),voxelResponse(stim1+nStimuli,:),'k.');
xlabel('Response to first presentation');
ylabel('Response to second presentation');
axis square
 
subplot(1,3,3);
plot(voxelResponse(stim1,:),voxelResponse(stim2,:),'k.');
xlabel(sprintf('Response to stimulus: %i deg',stimuli(stim1)));
ylabel(sprintf('Response to stimulus: %i',stimuli(stim2)));
axis square
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Make encoding model
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% make channel basis functions
nChannels = 8;
exponent = 7;
orientations = 0:179;
prefOrientation = 0:180/nChannels:179;
% loop over each channel
for iChannel = 1:nChannels
  % get sinusoid. Note the 2 here which makes it so that our 0 - 180 orientation
  % space gets mapped to all 360 degrees
  thisChannelBasis =  cos(2*pi*(orientations-prefOrientation(iChannel))/180);
  % rectify
  thisChannelBasis(thisChannelBasis<0) = 0;
  % raise to exponent
  thisChannelBasis = thisChannelBasis.^exponent;
  % keep in matrix
  channelBasis(:,iChannel) = thisChannelBasis;
end
 
% plot channel basis functions
figure;
plot(orientations,channelBasis);
xlabel('Preferred orientation (deg)');
ylabel('Ideal channel response (normalized to 1)');
 
% compute the channelResponse for each trial
for iTrial = 1:nTrials
  channelResponse(iTrial,:) = channelBasis(stimuli(iTrial),:);
end
 
% compute estimated weights
estimatedWeights =  pinv(channelResponse) * voxelResponse;
 
% compute model prediction
modelPrediction = channelResponse * estimatedWeights;
% compute residual
residualResponse = voxelResponse-modelPrediction;
% compute r2
r2 = 1-var(residualResponse(:))/var(voxelResponse(:))
 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Inverted encoding model
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% split half into train and test
firstHalf = 1:round(nTrials/2);
secondHalf = round(nTrials/2)+1:nTrials;
trainVoxelResponse = voxelResponse(firstHalf,:);
testVoxelResponse = voxelResponse(secondHalf:end,:);
% compute weights on train data
estimatedWeights = pinv(channelResponse(firstHalf,:))*trainVoxelResponse;
 
% compute channel response from textVoxelResponses
estimatedChannelResponse = testVoxelResponse * pinv(estimatedWeights);
 
% plot channel responses
figure;colors = hsv(nStimuli);
for iStimuli = 1:nStimuli
  plot(prefOrientation,mean(estimatedChannelResponse(iStimuli:nStimuli:end,:),1),'-','Color',colors(iStimuli,:));
  hold on
end
xlabel('Channel orientation preference (deg)');
ylabel('Estimated channel response (percentile of max)');
title(sprintf('r2 = %0.4f',r2));
 
% Compute voxel response without noise
nTrials = nStimuli * nRepeats;
for iTrial = 1:nTrials
  % get the neural response to this stimulus, by indexing the correct column of the neuralResponse matrix
  thisNeuralResponse = neuralResponse(:,stimuli(iTrial));
  % multiply this by the neuronToVoxelWeights to get the voxel response on this trial. Note that you need
  % to get the matrix dimensions right, so transpose is needed on thisNeuralResponse
  voxelResponseNoisy(iTrial,:) = thisNeuralResponse' * neuronToVoxelWeights;
end
 
% add noise
noiseStandardDeviation = 0.5;
% normalize response 
voxelResponseNoisy = voxelResponseNoisy / mean(voxelResponseNoisy(:));
% add gaussian noise
voxelResponseNoisy = voxelResponseNoisy + noiseStandardDeviation*randn(size(voxelResponseNoisy));
 
% split into train and tes
trainVoxelResponseNoisy = voxelResponseNoisy(firstHalf,:);
testVoxelResponseNoisy = voxelResponseNoisy(secondHalf:end,:);
 
% compute weights on train data
estimatedWeights = pinv(channelResponse(firstHalf,:))*trainVoxelResponseNoisy;
 
% compute model prediction on test data
modelPrediction = channelResponse(secondHalf,:) * estimatedWeights;
% compute residual
residualResponse = testVoxelResponseNoisy-modelPrediction;
% compute r2
r2 = 1-var(residualResponse(:))/var(testVoxelResponseNoisy(:))
 
% invert model and compute channel response
estimatedChannelResponse = testVoxelResponseNoisy * pinv(estimatedWeights);
 
% plot estimated channel profiles
figure;colors = hsv(nStimuli);
for iStimuli = 1:nStimuli
  plot(prefOrientation,mean(estimatedChannelResponse(iStimuli:nStimuli:end,:),1),'-','Color',colors(iStimuli,:));
  hold on
end
xlabel('Channel orientation preference (deg)');
ylabel('Estimated channel response (percentile of max)');
title(sprintf('r2 = %0.4f',r2));
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% stimulus likelihood function
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% split half into train and test
firstHalf = 1:round(nTrials/2);
secondHalf = round(nTrials/2)+1:nTrials;
trainVoxelResponse = voxelResponse(firstHalf,:);
testVoxelResponse = voxelResponse(secondHalf:end,:);
 
% compute weights on train data
estimatedWeights = pinv(channelResponse(firstHalf,:))*trainVoxelResponse;
 
% compute model prediction on test data
modelPrediction = channelResponse(secondHalf,:) * estimatedWeights;
% compute residual
residualResponse = testVoxelResponseNoisy-modelPrediction;
% compute residual variance, note that this is a scalar
residualVariance = var(residualResponse(:));
 
% make this into a covariance matrix in which the diagonal contains the variance for each voxel
% and off diagonals (in this case all 0) contain covariance between voxels
modelCovar = eye(nVoxels)*residualVariance;
 
% cycle over each trial
nTestTrials = size(testVoxelResponse,1);
for iTrial = 1:nTestTrials
  % now cycle over all possible orientation
  for iOrientation = 1:179
    % compute the mean voxel response predicted by the channel encoding model
    predictedResponse = channelBasis(iOrientation,:)*estimatedWeights;
    % now use that mean response and the model covariance to estimate the probability
    % of seeing this orientation given the response on this trial
    likelihood(iTrial,iOrientation) = mvnpdf(testVoxelResponse(iTrial,:),predictedResponse,modelCovar);
  end
end
 
figure
for iStimuli = 1:nStimuli
  plot(1:179,mean(likelihood(iStimuli:nStimuli:end,:),1),'-','Color',colors(iStimuli,:));
  hold on
end
xlabel('stimulus orientation (deg)');
ylabel('probability given trial response');
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Inverted Encoding model with different channel basis functions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% reweight the channels
channelReweighting = [0 0.8 0.4 0 0 0 0.4 0.8]';
% make into a full matrix xform to transform the original channels
for iChannel = 1:nChannels
  xform(iChannel,:) = circshift(channelReweighting,iChannel-1);
end
% and get new bimodal channels
bimodalChannelBasis = channelBasis * xform;
 
% display a figure with one of the channels
figure
plot(orientations,bimodalChannelBasis(:,5));
xlabel('orientation (deg)');
ylabel('Channel response (normalized to 1)');
 
% compute the channelResponse for each trial
for iTrial = 1:nTrials
  channelResponse(iTrial,:) = bimodalChannelBasis(stimuli(iTrial),:);
end
 
% compute estimated weights
estimatedWeights =  pinv(channelResponse) * voxelResponse;
 
% compute model prediction
modelPrediction = channelResponse * estimatedWeights;
% compute residual
residualResponse = voxelResponse-modelPrediction;
% compute r2
r2 = 1-var(residualResponse(:))/var(voxelResponse(:))
 
% compute estimated channel response profiles
estimatedChannelResponse = testVoxelResponse * pinv(estimatedWeights);
 
% and plot one of the channels averaged across all trials
figure;
plot(prefOrientation,mean(estimatedChannelResponse(5:nStimuli:end,:),1));
xlabel('Channel preferred orientation (deg)');
ylabel('Estimated channel response (percentile of full)');
title(sprintf('r2 = %0.4f',r2));
 
% compute weights on train data
estimatedWeights = pinv(channelResponse(firstHalf,:))*trainVoxelResponse;
 
% compute model prediction on test data
modelPrediction = channelResponse(secondHalf,:) * estimatedWeights;
% compute residual
residualResponse = testVoxelResponseNoisy-modelPrediction;
% compute residual variance, note that this is a scalar
residualVariance = var(residualResponse(:));
 
% make this into a covariance matrix in which the diagonal contains the variance for each voxel
% and off diagonals (in this case all 0) contain covariance between voxels
modelCovar = eye(nVoxels)*residualVariance;
 
% cycle over each trial
nTestTrials = size(testVoxelResponse,1);
for iTrial = 1:nTestTrials
  % now cycle over all possible orientation
  for iOrientation = 1:179
    % compute the mean voxel response predicted by the channel encoding model
    predictedResponse = bimodalChannelBasis(iOrientation,:)*estimatedWeights;
    % now use that mean response and the model covariance to estimate the probability
    % of seeing this orientation given the response on this trial
    likelihood(iTrial,iOrientation) = mvnpdf(testVoxelResponse(iTrial,:),predictedResponse,modelCovar);
  end
end
 
% now plot the likelihood function averaged over repeats
figure
for iStimuli = 1:nStimuli
  plot(1:179,mean(likelihood(iStimuli:nStimuli:end,:),1));
  hold on
end
xlabel('stimulus orientation (deg)');
ylabel('probability given trial response');