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

[pointprocesses] cleaned up code

Browse Source
This commit is contained in:
Jan Benda
2021-01-19 13:51:22 +01:00
parent 77b6666347
commit 25dd6ac2c7
12 changed files with 0 additions and 313 deletions
Download Patch File Download Diff File Expand all files Collapse all files

10
pointprocesses/code/colorednoisepdf.m
View File

@@ -1,10 +0,0 @@
function pcn = colorednoisepdf( x, misi, epsilon, tau )
% returns the ISI pdf for PIF with colored noise drive
% x: the input ISI
% misis: the mean isi
% epsilon: a parameter
% tau: the correlation time of the noise
gamma1 = x/tau+exp(-x/tau)-1.0;
gamma2 = 1.0-exp(-x/tau);
pcn=exp(-(x-misi).^2./(4.0*epsilon*tau.^2.*gamma1)).*(((misi-x).*gamma2+2*gamma1*tau).^2./(2*gamma1*tau^2)-epsilon*(gamma2.^2-2*gamma1.*exp(-x/tau))) ./ (2*tau*sqrt(4*pi*epsilon*gamma1.^3));
end

24
pointprocesses/code/colorednoisepdfisifit.m
View File

@@ -1,24 +0,0 @@
% misi = 0.02;
% epsilon = 1.0;
% tau = 0.1;
x=0:0.002:0.1;
% pcn = colorednoisepdf( x, misi, epsilon, tau )+10.0*randn( size( x ) );
% plot( x, pcn );
spikes = lifouspikes( 10, 15, 50.0, 1.0, 1.0 );
isivec = isis( spikes );
misi = mean( isivec );
1.0/misi
isibins = 0:0.0005:0.1;
[ n, c ] = hist( isivec, isibins );
n = n / sum(n)/(isibins(2)-isibins(1));
bar( c, n );
beta0 = [ 1.0, 0.01 ];
b = nlinfit(c(1:end-2), n(1:end-2), @(b,x)(colorednoisepdf(x, misi, b(1), b(2))), beta0)
pcn = colorednoisepdf( x, misi, b(1), b(2) );
hold on
plot( x, pcn, 'r', 'LineWidth', 3 );
hold off

55
pointprocesses/code/counthistpoisson.m
View File

@@ -1,55 +0,0 @@
function [counts, bins] = counthist(spikes, w)
% computes count histogram and compare with Poisson distribution
%
% [counts, bins] = counthist(spikes, w)
%
% Arguments:
% spikes: a cell array of vectors of spike times in seconds
% w: observation window duration in seconds for computing the counts
%
% Returns:
% counts: the histogram of counts normalized to probabilities
% bins: the bin centers for the histogram
% collect spike counts:
tmax = spikes{1}(end);
n = [];
r = [];
for k = 1:length(spikes)
times = spikes{k};
% alternative 1: count the number of spikes in each window:
% for tk = 0:w:tmax-w
% nn = sum( ( times >= tk ) & ( times < tk+w ) );
% %nn = length( find( ( times >= tk ) & ( times < tk+w ) ) );
% n = [ n nn ];
% end
% alternative 2: use the hist function to do that!
tbins = 0.5*w:w:tmax-0.5*w;
nn = hist(times, tbins);
n = [ n nn ];
% the rate of the spikes:
rate = (length(times)-1)/(times(end) - times(1));
r = [ r rate ];
end
% histogram of spike counts:
maxn = max( n );
[counts, bins ] = hist( n, 0:1:maxn+10 );
% normalize to probabilities:
counts = counts / sum( counts );
% plot:
if nargout == 0
bar( bins, counts );
hold on;
% Poisson distribution:
rate = mean( r );
x = 0:1:maxn+10;
a = rate*w;
y = a.^x.*exp(-a)./factorial(x);
plot( x, y, 'r', 'LineWidth', 3 );
hold off;
xlabel( 'counts k' );
ylabel( 'P(k)' );
end
end

24
pointprocesses/code/isireturnmap.m
View File

@@ -1,24 +0,0 @@
function isireturnmap( isis, lag2 )
% plot return maps for lag 1 and lag lag2
clf;
subplot( 1, 2, 1 );
lag = 1;
scatter( 1000.0*isis(1:end-lag)', 1000.0*isis(1+lag:end)', 'b', 'filled', 'MarkerEdgeColor', 'white' );
xlabel( 'ISI T_i [ms]' );
ylabel( 'ISI T_{i+1} [ms]' );
maxisi = max( isis );
maxy = ceil(maxisi/10)*10.0;
xlim( [0 1.5*maxy ])
ylim( [0 maxy ])
subplot( 1, 2, 2 );
lag = lag2;
scatter( 1000.0*isis(1:end-lag)', 1000.0*isis(1+lag:end)', 'b', 'filled', 'MarkerEdgeColor', 'white' );
xlabel( 'ISI T_i [ms]' );
ylabel( 'ISI T_{i+2} [ms]' );
xlim( [0 1.5*maxy ])
ylim( [0 maxy ])
end

100
pointprocesses/code/plotspikestats.m
View File

@@ -1,100 +0,0 @@
%% load data:
clear all
% alternative 1:
% pro: no structs. contra: global unknown variables
load poisson.mat
whos
poissonspikes = spikes;
load pifou.mat;
pifouspikes = spikes;
load lifadapt.mat;
lifadaptspikes = spikes;
clear spikes;
% alternative 2:
% pro: clean code. contra: structs that we do not really know yet
clear all
x = load( 'poisson.mat' );
poissonspikes = x.spikes;
x = load( 'pifou.mat' );
pifouspikes = x.spikes;
x = load( 'lifadapt.mat' );
lifadaptspikes = x.spikes;
%% spike raster plots:
tmax = 1.0;
subplot(1, 3, 1);
spikeraster(poissonspikes, tmax);
title('Poisson');
subplot(1, 3, 2);
spikeraster(pifouspikes, tmax);
title('PIF OU');
subplot(1, 3, 3);
spikeraster(lifadaptspikes, tmax);
title('LIF adapt');
%% isi histograms:
maxisi = 300.0;
binwidth = 0.002;
subplot(1, 3, 1);
poissonisis = isis(poissonspikes);
isihist(poissonisis, binwidth);
xlim([0, maxisi])
title('Poisson');
subplot(1, 3, 2);
pifouisis = isis(pifouspikes);
isihist(pifouisis, binwidth);
xlim([0, maxisi])
title('PIF OU');
subplot(1, 3, 3);
lifadaptisis = isis(lifadaptspikes);
isihist(lifadaptisis, binwidth);
xlim([0, maxisi])
title('LIF adapt');
%% serial correlations:
maxlag = 10;
rrange = [-0.5, 1.05];
subplot(1, 3, 1);
isiserialcorr(poissonisis, maxlag);
ylim(rrange)
title('Poisson');
subplot(1, 3, 2);
isiserialcorr(pifouisis, maxlag);
ylim(rrange)
title('PIF OU');
subplot(1, 3, 3);
isiserialcorr(lifadaptisis, maxlag);
ylim(rrange)
title('LIF adapt');
%% spike counts:
w = 0.1;
cmax = 8;
pmax = 0.5;
subplot(1, 3, 1);
counthist(poissonspikes, w);
xlim([0 cmax])
set(gca, 'XTick', 0:2:cmax)
ylim([0 pmax])
title('Poisson');
subplot(1, 3, 2);
counthist(pifouspikes, w);
xlim([0 cmax])
set(gca, 'XTick', 0:2:cmax)
ylim([0 pmax])
title('PIF OU');
subplot(1, 3, 3);
counthist(lifadaptspikes, w);
xlim([0 cmax])
set(gca, 'XTick', 0:2:cmax)
ylim([0 pmax])
title('LIF adapt');
savefigpdf(gcf, 'counthist.pdf', 20, 7);

27
pointprocesses/code/poissonisih.m
View File

@@ -1,27 +0,0 @@
rate = 100.0;
trials = 50;
tmax = 100.0;
% generate spikes:
spikes = poissonspikes( trials, rate, tmax );
% interspike intervals:
isivec = isis( spikes );
% histogram
f = figure( 1 );
isihist( isivec );
hold on
% theoretical density:
xmax = 5.0/rate;
x = 0:0.0001:xmax;
y = rate*exp(-rate*x);
plot( 1000.0*x, y, 'r', 'LineWidth', 3 );
% plot details:
title( sprintf( 'Poisson spike trains, rate=%g Hz, nisi=%d', rate, length( isivec ) ) )
xlim( [ 0.0 1000.0*xmax ] )
ylim( [ 0.0 1.1*rate ] )
legend( 'data', 'poisson' )
hold off
% serial correlations:
f = figure( 2 );
isiserialcorr( isivec, 10 );

46
pointprocesses/code/poissonisistats.m
View File

@@ -1,46 +0,0 @@
rates = 1:1:100;
avisi = [];
sdisi = [];
cvisi = [];
for rate = rates
spikes = poissonspikes( 10, rate, 100.0 );
isivec = isis( spikes );
av = mean( isivec );
sd = std( isivec );
cv = sd/av;
avisi = [ avisi av ];
sdisi = [ sdisi sd ];
cvisi = [ cvisi cv ];
end
f = figure;
subplot( 1, 3, 1 );
scatter( rates, 1000.0*avisi, 'b', 'filled' );
hold on;
plot( rates, 1000.0./rates, 'r' );
hold off;
xlabel( 'Rate \lambda [Hz]' );
ylim( [ 0 1000 ] );
title( 'Mean ISI [ms]' );
legend( 'simulation', 'theory 1/\lambda' );
subplot( 1, 3, 2 );
scatter( rates, 1000.0*sdisi, 'b', 'filled' );
hold on;
plot( rates, 1000.0./rates, 'r' );
hold off;
xlabel( 'Rate \lambda [Hz]' );
ylim( [ 0 1000 ] )
title( 'Standard deviation ISI [ms]' );
legend( 'simulation', 'theory 1/\lambda' );
subplot( 1, 3, 3 );
scatter( rates, cvisi, 'b', 'filled' );
hold on;
plot( rates, ones( size( rates ) ), 'r' );
hold off;
xlabel( 'Rate \lambda [Hz]' );
ylim( [ 0 2 ] )
title( 'CV' );
legend( 'simulation', 'theory' );

14
pointprocesses/code/psth.m
View File

@@ -1,14 +0,0 @@
function p = psth(spikes, dt, tmax)
% plots a PSTH of the spikes with binwidth dt
t = 0.0:dt:tmax+dt;
p = zeros(1, length(t));
for k=1:length(spikes)
times = spikes{k};
[h, b] = hist(times, t);
p = p + h;
end
p = p/length(spikes)/dt;
t(end) = [];
p(end) = [];
plot(t, p);
end

12
pointprocesses/code/spikerate.m
View File

@@ -1,12 +0,0 @@
function r = spikerate(spikes, duration)
% returns the average spike rate of the spikes
% for the first duration seconds
% spikes: a cell array of vectors of spike times
rates = zeros(length(spikes),1);
for k = 1:length(spikes)
times = spikes{k};
rates(k) = sum(times<duration)/duration;
end
r = mean(rates);
end