% function decode_neurons.m
%
% NSCI 850 - Computational Approaches to Neuroscience
%
% optimal Bayesian decoding of movements from neural spikes
% data set: curtesy of Lee Miller and Jim Rebesco
%
% G. Blohm (2010)
% modified from K. Kording (2008)
%--------------------------------------------------------------------------

load datasetReal
% loads X: firing rate of 35 neurons
%       y: classification (left or right)

%% Calculate the mean and std to compute the likelihood p(N|Left) and p(N|Right) for each neuron
for i=1:size(X,2)
    muLeft(i) = mean(X(find(y==0),i));
    stdLeft(i) = std(X(find(y==0),i));
    muRight(i) = mean(X(find(y==1),i));
    stdRight(i) = std(X(find(y==1),i));
end

%% Classify
% For each trial and each neuron calculate the likelihood p(N|Left) and p(N|Right)
% then compute the posterior probabilities for each neuron
for i=1:size(X,1)       % each trial
    for j=1:size(X,2)   % each neuron
        pNgLeft(i,j)=1/(sqrt(2*pi)*stdLeft(j))*exp(-(X(i,j)-muLeft(j)).^2/(2*stdLeft(j)^2))*1000;  % the multiplication with a fixed constant has no effect ... done for numerical stability
        pNgRight(i,j)=1/(sqrt(2*pi)*stdRight(j))*exp(-(X(i,j)-muRight(j)).^2/(2*stdRight(j)^2))*1000; 
        pLeftgN(i,j)=pNgLeft(i,j)./(pNgLeft(i,j)+pNgRight(i,j)); % the priors are hypothesized to be uniform here and can thus be neglected
        pRightgN(i,j)=pNgRight(i,j)./(pNgLeft(i,j)+pNgRight(i,j));
    end
end

% Combine neurons (works only if observations are statistically independent): see lecture slide 19
for i=1:size(X,1)
    pLeft(i)=prod(pLeftgN(i,:));
    pRight(i)=prod(pRightgN(i,:));
end
pLeft = pLeft./(pLeft+pRight);
pRight = pRight./(pLeft+pRight);

% check how good the classification is
percentCorrect=sum((pLeft<0.5)==y')/length(y)

%% plot results
figure % neural activity
X1=X(find(y==0),:);
X2=X(find(y==1),:);
imagesc([[X1; X2] 100*y]); % last column indicates left or right
colorbar
axis([1 36 1 205])
xlabel('# neuron')
ylabel('# trial')
title('Average firing rate')

