73 lines
1.8 KiB
Matlab
73 lines
1.8 KiB
Matlab
%% check Pig implementation:
|
|
mu = 0.01;
|
|
D = 1.0;
|
|
dx=0.0001;
|
|
x = dx:dx:10.0;
|
|
|
|
for D = [0.1 1 10 100]
|
|
for mu = [0.01 0.03 0.1 .03]
|
|
%pig = exp(-(x-mu).^2/4.0/D./x./mu/mu)./sqrt(4*pi*D*x.^3.0);
|
|
pig = invgauss(x,mu,D);
|
|
P = sum(pig)*dx;
|
|
mig = sum(x.*pig)*dx;
|
|
qig = sum(x.*x.*pig)*dx;
|
|
vig = qig - mig*mig;
|
|
dig = 0.5*vig/mig^3;
|
|
fprintf('Integral=%.3f\n', P);
|
|
fprintf('Average=%.3f = %.3f\n', mig, mu);
|
|
%fprintf('Variance=%.3f\n', vig);
|
|
fprintf('Diffusion=%.3f = %.3f\n\n', dig, D);
|
|
|
|
plot(x, pig)
|
|
hold on
|
|
plot([mu mu], [0 max(pig)], 'k')
|
|
hold off
|
|
xmax = 3.0*mu;
|
|
if xmax < 0.1
|
|
xmax = 0.1;
|
|
end
|
|
xlim([0 xmax])
|
|
text(0.6, 0.8, sprintf('D=%.3g', D), 'units', 'normalized')
|
|
pause( 1.0 )
|
|
end
|
|
end
|
|
|
|
|
|
%% fit to data:
|
|
trials = 1;
|
|
tmax = 100.0;
|
|
input = 10.0; % the input I
|
|
Dnoise = 1.0; % noise strength
|
|
outau = 0.1; % correlation time of the noise in seconds
|
|
%for input = [9.5 10.1 11.7 15.9]
|
|
% spikes = lifouspikes(trials, input, tmax, Dnoise, outau);
|
|
for input = [1.0 2.0 5.0 10.0]
|
|
spikes = pifouspikes(trials, input, tmax, Dnoise, outau);
|
|
isis = diff(spikes{1});
|
|
mu = mean(isis);
|
|
fprintf('mean firing rate = %.0f Hz\n', 1.0/mu)
|
|
v = var(isis);
|
|
D = 0.5*v/mu^3;
|
|
[h,b] = hist(isis, 0:0.001:1);
|
|
h = h/sum(h)/(b(2)-b(1));
|
|
bar(b, h);
|
|
hold on;
|
|
dx=0.0001;
|
|
x = dx:dx:max(isis);
|
|
pig = exp(-(x-mu).^2/4.0/D./x./mu/mu)./sqrt(4*pi*D*x.^3.0);
|
|
plot(x, pig, 'r', 'linewidth', 2);
|
|
% mle fit:
|
|
phat = mle(isis, 'pdf', @invgauss, 'start', [0.1 0.1] );
|
|
mu = phat(1);
|
|
D = phat(2);
|
|
pig = exp(-(x-mu).^2/4.0/D./x./mu/mu)./sqrt(4*pi*D*x.^3.0);
|
|
plot(x, pig, 'c', 'linewidth', 2);
|
|
xlim([0 4.0*mu])
|
|
hold off;
|
|
pause( 4.0 )
|
|
end
|
|
|
|
|
|
|
|
|