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

Handling of underscores in code filenames. Added code to pointprocesses.

Browse Source
This commit is contained in:
Jan Benda
2015-11-21 13:20:38 +01:00
parent 2fa2b4eab8
commit 788d6aa8ab
10 changed files with 141 additions and 61 deletions
Download Patch File Download Diff File Expand all files Collapse all files

9
pointprocesses/code/counthist.m
View File

@@ -1,9 +1,13 @@
function [counts, bins] = counthist(spikes, w)
% computes count histogram and compare them with Poisson distribution
% 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
@@ -27,11 +31,14 @@ function [counts, bins] = counthist(spikes, w)
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;

29
pointprocesses/code/isi_hist.m Normal file
View File

@@ -0,0 +1,29 @@
function [pdf, centers] = isi_hist(isis, binwidth)
% Compute normalized histogram of interspike intervals.
%
% [pdf, centers] = isi_hist(isis, binwidth)
%
% Arguments:
% isis: vector of interspike intervals in seconds
% binwidth: optional width in seconds to be used for the isi bins
%
% Returns:
% pdf: vector with probability density of interspike intervals in Hz
% centers: vector with corresponding centers of interspikeintervalls in seconds
if nargin < 2
% compute good binwidth:
nperbin = 200; % average number of data points per bin
bins = length( isis )/nperbin; % number of bins
binwidth = max( isis )/bins;
if binwidth < 5e-4 % half a millisecond
binwidth = 5e-4;
end
end
bins = 0.5*binwidth:binwidth:max(isis);
% histogram data:
[ nelements, centers ] = hist(isis, bins);
% normalization (integral = 1):
pdf = nelements / sum(nelements) / binwidth;
end

34
pointprocesses/code/isihist.m
View File

@@ -1,34 +0,0 @@
function isihist( isis, binwidth )
% plot histogram of interspike intervals
%
% isihist(isis, binwidth)
% isis: vector of interspike intervals in seconds
% binwidth: optional width in seconds to be used for the isi bins
if nargin < 2
nperbin = 200; % average number of data points per bin
bins = length( isis )/nperbin; % number of bins
binwidth = max( isis )/bins;
if binwidth < 5e-4 % half a millisecond
binwidth = 5e-4;
end
end
bins = 0.5*binwidth:binwidth:max(isis);
% histogram:
[ nelements, centers ] = hist( isis, bins );
% normalization (integral = 1):
nelements = nelements / sum( nelements ) / binwidth;
% plot:
bar( 1000.0*centers, nelements );
xlabel( 'ISI [ms]' )
ylabel( 'p(ISI) [1/s]')
% annotation:
misi = mean( isis );
sdisi = std( isis );
disi = sdisi^2.0/2.0/misi^3;
text( 0.95, 0.8, sprintf( 'mean=%.1f ms', 1000.0*misi ), 'Units', 'normalized', 'HorizontalAlignment', 'right' )
text( 0.95, 0.7, sprintf( 'std=%.1f ms', 1000.0*sdisi ), 'Units', 'normalized', 'HorizontalAlignment', 'right' )
text( 0.95, 0.6, sprintf( 'CV=%.2f', sdisi/misi ), 'Units', 'normalized', 'HorizontalAlignment', 'right' )
%text( 0.5, 0.3, sprintf( 'D=%.1f Hz', disi ), 'Units', 'normalized' )
end

10
pointprocesses/code/isiserialcorr.m
View File

@@ -1,9 +1,13 @@
function isicorr = isiserialcorr( isis, maxlag )
function isicorr = isiserialcorr(isis, maxlag)
% serial correlation of interspike intervals
%
% isicorr = isiserialcorr( isis, maxlag )
% isicorr = isiserialcorr(isis, maxlag)
%
% Arguments:
% isis: vector of interspike intervals in seconds
% maxlag: the maximum lag in seconds
%
% Returns:
% isicorr: vector with the serial correlations for lag 0 to maxlag
lags = 0:maxlag;
@@ -11,7 +15,7 @@ function isicorr = isiserialcorr( isis, maxlag )
for k = 1:length(lags)
lag = lags(k);
if length( isis ) > lag+10 % ensure "enough" data
% DANGER: the arguments to corr must be column vectors!
% NOTE: the arguments to corr must be column vectors!
% We insure this in the isis() function that generats the isis.
isicorr(k) = corr( isis(1:end-lag), isis(lag+1:end) );
end

28
pointprocesses/code/plot_isi_hist.m Normal file
View File

@@ -0,0 +1,28 @@
function plot_isi_hist(isis, binwidth)
% Plot and annotate histogram of interspike intervals.
%
% isihist(isis, binwidth)
%
% Arguments:
% isis: vector of interspike intervals in seconds
% binwidth: optional width in seconds to be used for the isi bins
if nargin < 2
[pdf, centers] = isi_hist(isis);
else
[pdf, centers] = isi_hist(isis, binwidth);
end
bar(1000.0*centers, nelements); % milliseconds on x-axis
xlabel('ISI [ms]')
ylabel('p(ISI) [1/s]')
% annotation:
misi = mean(isis);
sdisi = std(isis);
disi = sdisi^2.0/2.0/misi^3;
text(0.95, 0.8, sprintf('mean=%.1f ms', 1000.0*misi), 'Units', 'normalized', 'HorizontalAlignment', 'right')
text(0.95, 0.7, sprintf('std=%.1f ms', 1000.0*sdisi), 'Units', 'normalized', 'HorizontalAlignment', 'right')
text(0.95, 0.6, sprintf('CV=%.2f', sdisi/misi), 'Units', 'normalized', 'HorizontalAlignment', 'right')
%text(0.5, 0.3, sprintf('D=%.1f Hz', disi), 'Units', 'normalized')
end

13
pointprocesses/code/poissonspikes.m
View File

@@ -1,10 +1,14 @@
function spikes = poissonspikes( trials, rate, tmax )
% Generate spike times of a homogeneous poisson process
function spikes = poissonspikes(trials, rate, tmax)
% Generate spike times of a homogeneous poisson process.
%
% spikes = poissonspikes( trials, rate, tmax )
% spikes = poissonspikes(trials, rate, tmax)
%
% Arguments:
% trials: number of trials that should be generated
% rate: the rate of the Poisson process in Hertz
% tmax: the duration of each trial in seconds
%
% Returns:
% spikes: a cell array of vectors of spike times in seconds
dt = 3.33e-5;
@@ -16,7 +20,8 @@ function spikes = poissonspikes( trials, rate, tmax )
end
spikes = cell(trials, 1);
for k=1:trials
x = rand(round(tmax/dt), 1); % uniform random numbers for each bin
% uniform random numbers for each bin:
x = rand(round(tmax/dt), 1);
spikes{k} = find(x < p) * dt;
end
end