Merge branch 'master' of whale.am28.uni-tuebingen.de:scientificComputing

projects/project_mutualinfo/mutualinfo.tex

@@ -19,26 +19,34 @@
   object was reported by the subject. 
   
   \begin{parts}
-    \part Plot the data appropriately. 
+    \part Plot the data appropriately.
+
     \part Compute a 2-d histogram that shows how often different
-    combinations of reported and presented came up. 
+    combinations of reported and presented came up.
+
     \part Normalize the histogram such that it sums to one (i.e. make
     it a probability distribution $P(x,y)$ where $x$ is the presented
     object and $y$ is the reported object). Compute the probability
     distributions $P(x)$ and $P(y)$ in the same way.
+
     \part Use that probability distribution to compute the mutual
-    information $$I[x:y] = \sum_{x\in\{1,2\}}\sum_{y\in\{1,2\}} P(x,y)
+    information 
-    \log_2\frac{P(x,y)}{P(x)P(y)}$$ that the answers provide about the
+    \[ I[x:y] = \sum_{x\in\{1,2\}}\sum_{y\in\{1,2\}} P(x,y)
-    actually presented object. 
+    \log_2\frac{P(x,y)}{P(x)P(y)}\] 
+    that the answers provide about the actually presented object.
 
     The mutual information is a measure from information theory that is
     used in neuroscience to quantify, for example, how much information
     a spike train carries about a sensory stimulus.
+
     \part What is the maximally achievable mutual information (try to
     find out by generating your own dataset which naturally should
     yield maximal information)?
-    \part Use bootstrapping to compute the $95\%$ confidence interval
+
-    for the mutual information estimate in that dataset. 
+    \part Use bootstrapping (permutation test) to compute the $95\%$
+    confidence interval for the mutual information estimate in the
+    dataset from {\tt decisions.mat}.
+
   \end{parts}
 
 \end{questions}

projects/project_onset_fi/solution/boltzmannFit.m

@@ -0,0 +1,83 @@
+function [best_params] = boltzmannFit(contrasts, onset_fi)
+% use a gradient descent approach to fit the fi curve owith a Boltzmann
+% sigmoidal function
+%
+% Arguments:
+% constrasts: the x-values (must be sorted ascendingly)
+% onset_fi: the measured neuronal responses
+%
+% Returns the parameters leading to the best fit.
+
+% we will make some educated guesses for the start parameters; minimum:
+% minimum fi response; saturation: maximum fi response; inflection: zero
+% contrast; slope: response range/contrast range
+start_params = [min(onset_fi), max(onset_fi), 0, ...
+    (max(onset_fi) - min(onset_fi))/(max(contrasts) - min(contrasts))];
+
+best_params = gradientDescent(start_params, contrasts, onset_fi);
+
+
+%% subfunctions that could be separate m-files, but are only used in this context
+
+    function bp = gradientDescent(p, x, y)
+    % performs the gradient descent. Input arguments: 
+    % p: the list of starting parameters for the search
+    % x: the x-values
+    % y: the measured y-values
+    %
+    % returns the parameters leading to the best fit
+        
+       delta = 0.001;
+       gradient = getGradient(p, x, y, delta);
+       eps = [0.001 0.001, 0.00001, 0.00001]; % key point; use smaller changes for the more delicate parameters
+       while norm(gradient) > 0.5
+           p = p - eps .* gradient;
+           gradient = getGradient(p, x, y, delta);
+       end  
+       bp = p;
+    end
+
+    function g = getGradient(p, x, y, delta)
+    % calculate the gradient in the error surface for a small step
+    % along all parameters
+    % p: the parameter list
+    % x: the x-values
+    % y: the measured y-values
+    % delta: the step used or the gradient
+    %
+    % returns the gradient
+        g = zeros(size(p));
+        
+        for i = 1:length(p)
+            np = p;
+            np(i) = np(i) + delta;
+            g(i) = (costFunction(np, x, y) - costFunction(p, x, y)) / delta;
+        end
+    end
+
+    function cost = costFunction(p, x, y)
+        % cost function that relates the model to the data. Applies a
+        % Boltzman model and calculates the error in a mean square error
+        % sense.
+        % p: the set of parameters used for the model
+        % x: the x-values
+        % y: the measured y-values
+        %
+        % returns the costs related to that parameterization
+        
+        y_est = boltzmannModel(p, x);
+        cost = msqError(y, y_est);
+    end
+
+    function e = msqError(y, y_est)
+    % calculates the deviation between y and y_est in the mean square
+    % error sense
+    % y and y_est must have the same size
+    %
+    % returns the error
+        e = mean((y - y_est).^2);
+    end
+
+    
+
+end

projects/project_onset_fi/solution/boltzmannModel.m

@@ -0,0 +1,10 @@
+function y = boltzmannModel(p, x)
+% computes the Boltmannfunction given the parameters given in p
+% the following order is assumed
+% p(1) = minimum y-value
+% p(2) = maximum y-value, i.e. the saturation
+% p(3) = the position of the inflection point
+% p(4) = the slope at the inflection
+
+y = ((p(1) - p(2))) ./ (1 + exp((x-p(3))/p(4))) + p(2);
+end

projects/project_onset_fi/solution/convolutionRate.m

@@ -0,0 +1,27 @@
+function [time, rate] = convolution_rate(spikes, sigma, dt, t_max)
+% PSTH computed with convolution method. 
+%
+% [time, rate] = convolution_rate(spikes, sigma, dt, t_max)
+%
+% Arguments: 
+%   spikes: a vector containing the spike times.
+%   sigma : the standard deviation of the Gaussian kernel in seconds.
+%   dt    : the temporal resolution in seconds.
+%   t_max : the trial duration in seconds.
+%
+% Returns:    two vectors containing the time and the rate.
+
+time = 0:dt:t_max - dt;
+rate = zeros(size(time));
+spike_indices = round(spikes / dt);
+rate(spike_indices) = 1;
+kernel = gaussKernel(sigma, dt);
+
+rate = conv(rate, kernel, 'same');
+end
+
+
+function y = gaussKernel(s, step)
+    x = -4 * s:step:4 * s;
+    y = exp(-0.5 .* (x ./ s) .^ 2) ./ sqrt(2 * pi) / s;
+end

projects/project_onset_fi/solution/getFICurve.m

@@ -0,0 +1,12 @@
+function [fi, fi_std, sort_contrasts] = getFICurve(firing_rates, time, duration, contrasts)
+
+[sort_contrasts, sort_idx] = sort(contrasts);
+fi = zeros(length(sort_contrasts), 1);
+fi_std = zeros(length(sort_contrasts), 1); 
+
+for i = 1:length(sort_contrasts)
+    responses = firing_rates{sort_idx(i)};
+    onset_responses = mean(responses((time > 0) & (time <= duration), :),1);
+    fi(i) = mean(onset_responses);
+    fi_std(i) = std(onset_responses); 
+end

projects/project_onset_fi/solution/getFiringRates.m

@@ -0,0 +1,10 @@
+function [time, rates] = getFiringRates(all_spikes, min_t, max_t, dt)
+
+
+time = min_t:dt:max_t-dt;
+rates = zeros(length(time), length(all_spikes));
+for i = 1:length(all_spikes)
+    spike_times = all_spikes{i}/1000 - min_t;
+    spike_times(spike_times == 0.0) = dt;
+    [t, rates(:, i)] = convolutionRate(spike_times, 0.01, dt, max_t-min_t);
+end

projects/project_onset_fi/solution/main.m

@@ -0,0 +1,36 @@
+%% solution for project onset-FI
+
+clear all 
+close all
+
+dt = 1./20000;
+t_min = -0.2;
+t_max = 1.825;
+data_folder = strcat('..', filesep, 'data');
+files = dir(strcat(data_folder, filesep, '*.mat'));
+
+contrasts = zeros(length(files), 1);
+rates = {};
+for i = 1:length(files)
+    data = load(strcat(data_folder, filesep, files(i).name));
+    contrasts(i) = data.contrast;
+    spikes_times = data.spike_times;
+    [t, rates{i}] = getFiringRates(spikes_times, t_min, t_max, dt);
+end
+
+%% create a plot of the average response for the different contrasts
+plotAverageResponse(rates, t, contrasts, 'averageResponse')
+
+%% extract the onset and the steady-state firing rate
+[fi, fi_std, contr] = getFICurve(rates, t, 0.015, contrasts);
+
+%% plot FI curve 
+plotFICurve(fi, fi_std, contr, 'ficurve')
+
+%% fit FICurve with Boltzmann function
+best_params = boltzmannFit(contr, fi);
+
+% plot the fi curve with the fit
+plotFICurveFit(fi, fi_std, contr, best_params, 'ficurvefit')
+
+

projects/project_onset_fi/solution/plotAverageResponse.m

@@ -0,0 +1,47 @@
+function [] = plotAverageResponse(rates, time, contrasts, filename)
+[sort_contrasts, sort_idx] = sort(contrasts);
+
+figure
+set(gcf, 'paperunits', 'centimeters', 'papersize', [20 20], ...
+    'paperposition', [0.0 0.0 20, 20], 'color', 'white')
+axes = cell(size(sort_contrasts));
+for i = 1:length(sort_contrasts)
+    subplot(ceil(length(sort_contrasts)/2), 2, i);
+    responses = rates{sort_idx(i)};
+    avg_response = mean(responses, 2)';
+    std_response = std(responses, [], 2)';
+    f = fill([time fliplr(time)], [avg_response + std_response fliplr(avg_response - std_response)], ...
+            'b', 'EdgeColor','none', 'displayname', 'std');
+    alpha(f, 0.25)
+    hold on
+    plot(time, avg_response, 'b', 'displayname', 'average response')
+    text(1.2, 500, sprintf('contrast: %i %', sort_contrasts(i)), ...
+        'fontsize', 9, 'FontWeight', 'bold')
+    xlim([time(1) time(end)])
+    ylim([0, 600])
+    set(gca, 'XTick', time(1):0.2:time(end), 'YTick', 0:200:600)
+    if mod(i,2) == 1
+        ylabel('firing rate [Hz]', 'fontsize', 11);
+    else
+        set(gca, 'yticklabels', []);
+    end
+    if i > length(sort_contrasts) - 2 
+        xlabel('time [s]', 'fontsize', 11)
+    else
+        set(gca, 'xticklabels', []);
+    end
+   
+    box('off')
+    axes{i} = gca();
+    if i == 1
+        l = legend('show');
+        l.FontSize = 9;
+        l.Box = 'off';
+        l.Orientation = 'horizontal';
+        l.Position = [0.25, 0.96, .5, 0.01];
+    end
+end
+
+axes{1}.Position = [axes{3}.Position(1) axes{2}.Position(2) ... 
+    axes{2}.Position(3) axes{3}.Position(4)];
+saveas(gcf, filename, 'pdf')

projects/project_onset_fi/solution/plotFICurve.m

@@ -0,0 +1,16 @@
+function plotFICurve(fi, fi_std, contrasts, filename)
+
+figure()
+set(gcf, 'paperunits', 'centimeters', 'papersize', [10 10], ...
+    'paperposition', [0.0 0.0 10, 10], 'color', 'white')
+
+errorbar(contrasts, fi, fi_std, 'o', 'displayname', 'onset response')
+xlim([-20, 20])
+ylim([0, 600])
+xlabel('stimulus contrast [%]', 'fontsize', 11)
+ylabel('firing rate [Hz]', 'fontsize', 11)
+box('off')
+
+if ~isempty(filename) 
+    saveas(gcf, filename, 'pdf')
+end

projects/project_onset_fi/solution/plotFICurveFit.m

@@ -0,0 +1,8 @@
+function plotFICurveFit(fi, fi_error, contrasts, fit_parameter, filename)
+
+plotFICurve(fi, fi_error, contrasts, '')
+fit = boltzmannModel(fit_parameter, min(contrasts):0.1:max(contrasts));
+hold on 
+plot(min(contrasts):0.1:max(contrasts), fit, 'displayname', 'Boltzmann fit')
+legend('Location', 'northwest')
+saveas(gcf, filename, 'pdf')

projects/project_populationvector/solution/maxlikelihoodangle.m

@@ -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

projects/project_populationvector/solution/populationvector.m

@@ -1,24 +1,24 @@
 close all
-datapath = '../';
+datapath = '../data/';
-% datapath = '../code/';
 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:
 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'));
 phases = zeros(length(files), 1);
+gains = zeros(length(files), 1);
 figure()
 for j = 1:length(files)
 	file = files(j);
@@ -31,10 +31,10 @@ for j = 1:length(files)
 		rates(k) = r;
 	end
 	[mr, maxi] = max(rates);
-	p0 = [mr, angles(maxi)/180.0-0.5];
+	p0 = [mr, angles(maxi)];
 	%p = p0;
 	p = lsqcurvefit(cosine, p0, angles, rates');
-	phase = p(2)*180.0;
+	phase = p(2);
     if phase > 180.0
         phase = phase - 180.0;
     end
@@ -42,10 +42,12 @@ for j = 1:length(files)
         phase = phase + 180.0;
     end
     phases(j) = phase;
+    gains(j) = p(1);
 	subplot(2, 3, j);
 	plot(angles, rates, 'b');
 	hold on;
-	plot(angles, cosine(p, angles), 'r');
+    a = 0:0.1:180;
+	plot(a, cosine(p, a), 'r');
 	hold off;
 	xlim([0.0 180.0])
 	ylim([0.0 50.0])
@@ -57,12 +59,13 @@ end
 a = load(strcat(datapath, 'population04.mat'));
 spikes = a.spikes;
 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();
-subplot(1, 3, 1);
+subplot(2, 2, 1);
 angleestimates1 = zeros(size(spikes, 2), 1);
 angleestimates2 = zeros(size(spikes, 2), 1);
+angleestimates3 = zeros(size(spikes, 2), 1);
 [x, inx] = sort(phases);
 % loop over trials:
 for j = 1:size(spikes, 2)
@@ -76,35 +79,60 @@ for j = 1:size(spikes, 2)
     angleestimates1(j) = popvecangle(phases, rates);
     [m, i] = max(rates);
     angleestimates2(j) = phases(i);
+    angleestimates3(j) = maxlikelihoodangle(phases, gains, rates);
 end
 xlabel('preferred angle')
 ylabel('firing rate')
 hold off;
-subplot(1, 3, 2);
+subplot(2, 2, 2);
 hist(angleestimates1);
-xlabel('stimulus angle')
+xlabel('population vector angle')
-subplot(1, 3, 3);
+subplot(2, 2, 3);
 hist(angleestimates2);
-xlabel('stimulus angle')
+xlabel('max. rate angle')
+subplot(2, 2, 4);
+hist(angleestimates3);
+xlabel('max. likelihood angle')
 angle
 mean(angleestimates1)
 mean(angleestimates2)
+mean(angleestimates3)
 
 
 %% read out robustness:
 files = dir(strcat(datapath, 'population*.mat'));
 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)
 	file = files(i);
 	a = load(strcat(datapath, file.name));
     spikes = a.spikes;
     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);
     angleestimates2 = zeros(size(spikes, 2), 1);
+    angleestimates3 = zeros(size(spikes, 2), 1);
     for j = 1:size(spikes, 2)
         rates = zeros(size(spikes, 1), 1);
         for k = 1:size(spikes, 1)
@@ -114,19 +142,52 @@ for i = 1:length(files)
         angleestimates1(j) = popvecangle(phases, rates);
         [m, inx] = max(rates);
         angleestimates2(j) = phases(inx);
+        angleestimates3(j) = maxlikelihoodangle(phases, gains, rates);
     end
     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
+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();
-subplot(1, 2, 1);
+subplot(1, 3, 1);
-scatter(angles, e1m);
+hold on;
-xlabel('stimuluis angle')
+plot(x, x, 'k');
-ylabel('estimated angle')
+scatter(angles, e1sm);
-subplot(1, 2, 2);
+xlabel('stimulus angle')
-scatter(angles, e2m);
+ylabel('estimated angle (population vector)')
-xlabel('stimuluis angle')
+subplot(1, 3, 2);
-ylabel('estimated angle')
+hold on;
+plot(x, x, 'k');
+scatter(angles, e2sm);
+xlabel('stimulus angle')
+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)')