teaching/scientificComputing
Archived
3
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

added maximum likelihood solution

Browse Source
This commit is contained in:
Jan Benda 2018-01-30 14:37:28 +01:00
parent 2103216caf
commit 40dbc3cfb9
2 changed files with 116 additions and 31 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

23
projects/project_populationvector/solution/maxlikelihoodangle.m Normal file
View File

@ -0,0 +1,23 @@
function angle = maxlikelihoodangle(phases, gains, rates)
% maximum likelihood estimation of orientation
cosine = @(g,p,xdata)0.5*g.*(1.0+cos(2.0*pi*(xdata-p)/180.0));
angels = 0:1.0:180.0;
loglikelihoods = zeros(length(phases), length(angels));
for i=1:length(phases)
r = cosine(gains(i), phases(i), angels);
%%loglikelihoods(i, :) = exp(-0.5*((rates(i)-r)./(0.25*r)).^2.0)./sqrt(2.0*pi*(0.25*r).^2.0);
%loglikelihoods(i, :) = log(exp(-0.5*((rates(i)-r)./(0.25*r)).^2.0)./sqrt(2.0*pi*(0.25*r).^2.0));
loglikelihoods(i, :) = -0.5*((rates(i)-r)./(0.25*r)).^2.0 - 0.5*log(2.0*pi*(0.25*r).^2.0);
end
loglikelihood = sum(loglikelihoods, 1);
[m i] = max(loglikelihood);
angle = angels(i);
% plot(angels, loglikelihood);
% hold on;
% plot([angle angle], [-500 0], 'k')
% hold off;
% xlabel('angle');
% ylabel('likelihood');
% ylim([-500, 0]);
% pause( 0.2 );
end

124
projects/project_populationvector/solution/populationvector.m
View File

@ -1,23 +1,24 @@
close all close all
datapath = '../data/'; datapath = '../data/';
files = dir(strcat(datapath, 'unit*.mat')); files = dir(strcat(datapath, 'unit*.mat'));
for file = files' % for file = files'
a = load(strcat(datapath, file.name)); % a = load(strcat(datapath, file.name));
spikes = a.spikes; % spikes = a.spikes;
angles = a.angles; % angles = a.angles;
figure() % figure()
for k = 1:size(spikes, 1) % for k = 1:size(spikes, 1)
subplot(3, 4, k) % subplot(3, 4, k)
spikeraster(spikes(k,:), -0.2, 0.6); % spikeraster(spikes(k,:), -0.2, 0.6);
end % end
end % end
%% tuning curves: %% tuning curves:
close all close all
cosine = @(p,xdata)0.5*p(1).*(1.0+cos(2.0*pi*(xdata/180.0-p(2)))); cosine = @(p,xdata)0.5*p(1).*(1.0+cos(2.0*pi*(xdata-p(2))/180.0));
files = dir(strcat(datapath, 'unit*.mat')); files = dir(strcat(datapath, 'unit*.mat'));
phases = zeros(length(files), 1); phases = zeros(length(files), 1);
gains = zeros(length(files), 1);
figure() figure()
for j = 1:length(files) for j = 1:length(files)
file = files(j); file = files(j);
@ -30,10 +31,10 @@ for j = 1:length(files)
rates(k) = r; rates(k) = r;
end end
[mr, maxi] = max(rates); [mr, maxi] = max(rates);
p0 = [mr, angles(maxi)/180.0-0.5]; p0 = [mr, angles(maxi)];
%p = p0; %p = p0;
p = lsqcurvefit(cosine, p0, angles, rates'); p = lsqcurvefit(cosine, p0, angles, rates');
phase = p(2)*180.0; phase = p(2);
if phase > 180.0 if phase > 180.0
phase = phase - 180.0; phase = phase - 180.0;
end end
@ -41,10 +42,12 @@ for j = 1:length(files)
phase = phase + 180.0; phase = phase + 180.0;
end end
phases(j) = phase; phases(j) = phase;
gains(j) = p(1);
subplot(2, 3, j); subplot(2, 3, j);
plot(angles, rates, 'b'); plot(angles, rates, 'b');
hold on; hold on;
plot(angles, cosine(p, angles), 'r'); a = 0:0.1:180;
plot(a, cosine(p, a), 'r');
hold off; hold off;
xlim([0.0 180.0]) xlim([0.0 180.0])
ylim([0.0 50.0]) ylim([0.0 50.0])
@ -56,12 +59,13 @@ end
a = load(strcat(datapath, 'population04.mat')); a = load(strcat(datapath, 'population04.mat'));
spikes = a.spikes; spikes = a.spikes;
angle = a.angle; angle = a.angle;
unitphases = a.phases*180.0; % unitphases = a.phases*180.0;
unitphases(unitphases>180.0) = unitphases(unitphases>180.0) - 180.0; % unitphases(unitphases>180.0) = unitphases(unitphases>180.0) - 180.0;
figure(); figure();
subplot(1, 3, 1); subplot(2, 2, 1);
angleestimates1 = zeros(size(spikes, 2), 1); angleestimates1 = zeros(size(spikes, 2), 1);
angleestimates2 = zeros(size(spikes, 2), 1); angleestimates2 = zeros(size(spikes, 2), 1);
angleestimates3 = zeros(size(spikes, 2), 1);
[x, inx] = sort(phases); [x, inx] = sort(phases);
% loop over trials: % loop over trials:
for j = 1:size(spikes, 2) for j = 1:size(spikes, 2)
@ -75,35 +79,60 @@ for j = 1:size(spikes, 2)
angleestimates1(j) = popvecangle(phases, rates); angleestimates1(j) = popvecangle(phases, rates);
[m, i] = max(rates); [m, i] = max(rates);
angleestimates2(j) = phases(i); angleestimates2(j) = phases(i);
angleestimates3(j) = maxlikelihoodangle(phases, gains, rates);
end end
xlabel('preferred angle') xlabel('preferred angle')
ylabel('firing rate') ylabel('firing rate')
hold off; hold off;
subplot(1, 3, 2); subplot(2, 2, 2);
hist(angleestimates1); hist(angleestimates1);
xlabel('population vector angle') xlabel('population vector angle')
subplot(1, 3, 3); subplot(2, 2, 3);
hist(angleestimates2); hist(angleestimates2);
xlabel('max. rate angle') xlabel('max. rate angle')
subplot(2, 2, 4);
hist(angleestimates3);
xlabel('max. likelihood angle')
angle angle
mean(angleestimates1) mean(angleestimates1)
mean(angleestimates2) mean(angleestimates2)
mean(angleestimates3)
%% read out robustness: %% read out robustness:
files = dir(strcat(datapath, 'population*.mat')); files = dir(strcat(datapath, 'population*.mat'));
angles = zeros(length(files), 1); angles = zeros(length(files), 1);
e1m = zeros(length(files), 1); e1mm = zeros(length(files), 1);
e1s = zeros(length(files), 1); e2mm = zeros(length(files), 1);
e2m = zeros(length(files), 1); e3mm = zeros(length(files), 1);
e2s = zeros(length(files), 1); e1sm = zeros(length(files), 1);
e1ss = zeros(length(files), 1);
e2sm = zeros(length(files), 1);
e2ss = zeros(length(files), 1);
e3sm = zeros(length(files), 1);
e3ss = zeros(length(files), 1);
for i = 1:length(files) for i = 1:length(files)
file = files(i); file = files(i);
a = load(strcat(datapath, file.name)); a = load(strcat(datapath, file.name));
spikes = a.spikes; spikes = a.spikes;
angle = a.angle; angle = a.angle;
% multi trial estimates:
rates = zeros(size(spikes, 1), 1);
for k = 1:size(spikes, 1)
r = zeros(size(spikes, 2), 1);
for j = 1:size(spikes, 2)
r(j) = firingrate(spikes(k, j), 0.0, 0.2);
end
rates(k) = mean(r);
end
e1mm(i) = popvecangle(phases, rates);
[m, inx] = max(rates);
e2mm(i) = phases(inx);
e3mm(i) = maxlikelihoodangle(phases, gains, rates);
% single trial estimates:
angleestimates1 = zeros(size(spikes, 2), 1); angleestimates1 = zeros(size(spikes, 2), 1);
angleestimates2 = zeros(size(spikes, 2), 1); angleestimates2 = zeros(size(spikes, 2), 1);
angleestimates3 = zeros(size(spikes, 2), 1);
for j = 1:size(spikes, 2) for j = 1:size(spikes, 2)
rates = zeros(size(spikes, 1), 1); rates = zeros(size(spikes, 1), 1);
for k = 1:size(spikes, 1) for k = 1:size(spikes, 1)
@ -113,19 +142,52 @@ for i = 1:length(files)
angleestimates1(j) = popvecangle(phases, rates); angleestimates1(j) = popvecangle(phases, rates);
[m, inx] = max(rates); [m, inx] = max(rates);
angleestimates2(j) = phases(inx); angleestimates2(j) = phases(inx);
angleestimates3(j) = maxlikelihoodangle(phases, gains, rates);
end end
angles(i) = angle; angles(i) = angle;
e1m(i) = mean(angleestimates1); e1sm(i) = mean(angleestimates1);
e1s(i) = std(angleestimates1); e1ss(i) = std(angleestimates1);
e2m(i) = mean(angleestimates2); e2sm(i) = mean(angleestimates2);
e2s(i) = std(angleestimates2); e2ss(i) = std(angleestimates2);
e3sm(i) = mean(angleestimates3);
e3ss(i) = std(angleestimates3);
end end
x = 0:180.0;
figure();
subplot(1, 3, 1);
hold on;
plot(x, x, 'k');
scatter(angles, e1mm);
xlabel('stimulus angle')
ylabel('estimated angle (population vector)')
subplot(1, 3, 2);
hold on;
plot(x, x, 'k');
scatter(angles, e2mm);
xlabel('stimulus angle')
ylabel('estimated angle (maximum firing rate)')
subplot(1, 3, 3);
hold on;
plot(x, x, 'k');
scatter(angles, e3mm);
xlabel('stimulus angle')
ylabel('estimated angle (maximum likelihood)')
figure(); figure();
subplot(1, 2, 1); subplot(1, 3, 1);
scatter(angles, e1m); hold on;
plot(x, x, 'k');
scatter(angles, e1sm);
xlabel('stimulus angle') xlabel('stimulus angle')
ylabel('estimated angle (population vector)') ylabel('estimated angle (population vector)')
subplot(1, 2, 2); subplot(1, 3, 2);
scatter(angles, e2m); hold on;
plot(x, x, 'k');
scatter(angles, e2sm);
xlabel('stimulus angle') xlabel('stimulus angle')
ylabel('estimated angle (maximum firing rate)') ylabel('estimated angle (maximum firing rate)')
subplot(1, 3, 3);
hold on;
plot(x, x, 'k');
scatter(angles, e3sm);
xlabel('stimulus angle')
ylabel('estimated angle (maximum likelihood)')