% Test my matlab code
clear all
close all

% Set up model
model.pstim = 0.9; % probability of successfully eliciting AP in this terminal
model.offset = 63; % to be added to all EPSPs generated
model.n = 5; % number of release sites
model.q = 150; % contribution from each release site which releases a quantum
model.p = 0.4; % probability that each release site releases a quantum if stimulated
model.QVtype = 2; % 1=Type I, 2=flat
model.qsigma = 20; % different interpretations for type I or flat
model.noise = 40; % SD of Gaussian noise


% Generate simulated data:
ndata=1e5;
[simdata,pfailure] = Simulate(model,ndata);

% Bin simulata
[freq,v]=hist(simdata,50);
figure
bar(v,freq)
xlabel('voltage of EPSP')
ylabel('number of events')

binwidth=min(diff(v));

% Superimpose PDF:
vv = [min(v) : range(v)/1000 : max(v)];
PDF = BinomPDF(model,vv);
hold on
plot(vv,PDF*ndata*binwidth,'r');
% NB the *ndata*binwidth scales the PDF up so it matches height in the
% frequency bins


% Now add in a glitch rate.
model.glitchrate = 0.4;
model.glitchmin = -1200;
model.glitchmax = 2800;
% Say: each EPSP has a [model.glitchrate] chance of being replaced with a
% new value picked randomly from the interval [u1,u2]
newsimdata = simdata;
jreplace = find(rand(size(newsimdata))<model.glitchrate);
newsimdata(jreplace) = model.glitchmin + (model.glitchmax-model.glitchmin)*rand(size(jreplace)) + model.offset;
figure
[freq,v]=hist(newsimdata,80);
bar(v,freq)
binwidth=min(diff(v));
xlabel('voltage of EPSP')
ylabel('number of events')
vv = [min( [model.glitchmin min(v)]) : range(v)/10000 : max([model.glitchmax max(v)])];
PDF = BinomPDF(model,vv);
hold on
plot(vv,PDF*ndata*binwidth,'g');


