diff --git a/cont_ca_sampler.m b/cont_ca_sampler.m index adce0d7..58f1fb1 100644 --- a/cont_ca_sampler.m +++ b/cont_ca_sampler.m @@ -1,10 +1,9 @@ function SAMPLES = cont_ca_sampler(Y,params) -% add the capability of dealing with missing numbers +% MCMC continuous time sampler of spiketimes given a fluorescence trace Y -% Continuous time sampler +% Inputs: % Y data -% P intialization parameters (discrete time constant P.g required) % params parameters structure % params.g discrete time constant(s) (estimated if not provided) % params.sn initializer for noise (estimated if not provided) @@ -12,7 +11,7 @@ % params.c1 initializer for initial concentration (estimated if not provided) % params.Nsamples number of samples after burn in (default 500) % params.B number of burn in samples (default 200) -% params.marg flag for marginalized sampler (default 0) +% params.marg flag for marginalized sampler for baseline and initial concentration (default 0) % params.upd_gam flag for updating gamma (default 1) % params.gam_step number of samples after which gamma is updated (default 50) % params.std_move standard deviation of shifting kernel (default 3*Dt) @@ -45,6 +44,9 @@ % Author: Eftychios A. Pnevmatikakis and Josh Merel +% Reference: Pnevmatikakis et al., Bayesian spike inference from calcium +% imaging data, Asilomar 2013. + Y = Y(:); T = length(Y); isanY = ~isnan(Y); % Deal with possible missing entries diff --git a/plot_continuous_samples.m b/plot_continuous_samples.m index fbd3fe2..52a832a 100755 --- a/plot_continuous_samples.m +++ b/plot_continuous_samples.m @@ -1,5 +1,18 @@ function plot_continuous_samples(SAMPLES,Y) +% plot results of MCMC sampler +% The mean calcium sample, spike sampler raster plot and samples for +% amplitude, number of spikes, discrete time constants, noise variance, +% baseline and initial concentration are generated, together with their +% autocorrelation functions. If the marginalized flag was used, then the +% posterior pdfs of baseline and initial concentration are plotted. + +% Inputs: +% SAMPLES: structure with SAMPLES obtained from cont_ca_sampler.m +% Y: inpurt fluorescence trace + +% Author: Eftychios A. Pnevmatikakis, 2016, Simons Foundation + T = length(Y); N = length(SAMPLES.ns); show_gamma = 1; @@ -33,7 +46,6 @@ function plot_continuous_samples(SAMPLES,Y) end - if length(SAMPLES.Cb) == 2 marg = 1; % marginalized sampler else diff --git a/sampling_demo_ar2.m b/sampling_demo_ar2.m index feec0b9..43c885e 100644 --- a/sampling_demo_ar2.m +++ b/sampling_demo_ar2.m @@ -1,6 +1,8 @@ -clear all; +% creating a demo for the MCMC sampler + +clearvars; addpath utilities -dt = 1e-1; +dt = 1e-1; % Data is generated with provided dt but sampled with Dt = 1 T = 7000; ld = 0.05; % rate spikes per second @@ -10,22 +12,23 @@ hmax = tau_decay/(tau_decay+tau_rise)*(tau_rise/(tau_decay+tau_rise))^(tau_rise/tau_decay); [g,h1] = tau_c2d(tau_rise,tau_decay,dt); -b = hmax/4; -cin = [.2*b,.15*b]; -c = [cin,filter(h1,[1,-g],s(3:end),filtic(h1,[1,-g],cin))] + b; -c_true = c(round(1/dt):round(1/dt):round(T/dt)); -sg = hmax/4; % noise level -y = c_true + sg*randn(1,length(c_true)); +b = hmax/4; % baseline value +cin = [.2*b,.15*b]; % initial concentration +c = [cin,filter(h1,[1,-g],s(3:end),filtic(h1,[1,-g],cin))] + b; % true calcium at original resolution +c_true = c(round(1/dt):round(1/dt):round(T/dt)); % true calcium at sampled resolution +sg = hmax/4; % noise level +y = c_true + sg*randn(1,length(c_true)); % noisy observations figure;subplot(2,1,1); plot(dt:dt:T,c); hold all; scatter(1:T,y,'r*'); legend('True Calcium','Observed Values'); -%% constrained foopsi -[g2,h2] = tau_c2d(tau_rise,tau_decay,1); -P.f = 1; -P.g = g2; -P.sn = sg; -params.marg = 0; + +%% Run constrained deconvolution approach (constrained foopsi) + +[g2,h2] = tau_c2d(tau_rise,tau_decay,1); % generate discrete time constants + [ca_foopsi,cb,c1,~,~,spikes_foopsi] = constrained_foopsi(y,[],[],g2); + +% do some plotting spiketimes{1} = find(s)*dt; spikeRaster = samples_cell2mat(spiketimes,T,1); f = find(spikeRaster); @@ -44,7 +47,8 @@ title('Foopsi Spikes','FontWeight','bold','Fontsize',14); xlabel('Timestep','FontWeight','bold','Fontsize',16); legend('Foopsi Spikes','Ground Truth'); drawnow; -%% MCMC + +%% run continuous time MCMC sampler params.p = 2; params.g = g2; @@ -54,5 +58,7 @@ params.c1 = c1; params.sn = sg; params.marg = 0; -SAMPLES2 = cont_ca_sampler(y,params); %% MCMC -plot_continuous_samples(SAMPLES2,y(:)); \ No newline at end of file + +SAMPLES = cont_ca_sampler(y,params); %% MCMC +plot_continuous_samples(SAMPLES,y(:)); +M = plot_marginals(SAMPLES.ss,T,spikeRaster); \ No newline at end of file diff --git a/utilities/addSpike.m b/utilities/addSpike.m index 1f94044..db85eb4 100755 --- a/utilities/addSpike.m +++ b/utilities/addSpike.m @@ -1,12 +1,33 @@ -function [newSpikeTrain, newCalcium, newLL] = addSpike(oldSpikeTrain,oldCalcium,oldLL,filter,tau,obsCalcium,timeToAdd,indx,Dt,A) +function [newSpikeTrain, newCalcium, newLL] = addSpike(oldSpikeTrain,oldCalcium,oldLL,filters,tau,obsCalcium,timeToAdd,indx,Dt,A) + +% Add a given spike to the existing spike train. + +% Inputs: +% oldSpikeTrain: current spike train +% oldCalcium: current noiseless calcium trace +% oldLL: current value of the log-likelihood function +% filters: exponential rise and decay kernels for calcium transient +% tau: continuous time rise and decay time constants +% obsCalcium: observed fluorescence trace +% timetoAdd: time of the spike to be added +% indx: place where the new spike is added in the existing spike train vector +% Dt: time-bin width +% A: spike amplitude + +% Outputs: +% newSpikeTrain: new vector of spike times +% newCalcium: new noiseless calcium trace +% newLL: new value of the log-likelihood function + +% Author: Eftychios A. Pnevmatikakis and Josh Merel tau_h = tau(1); tau_d = tau(2); - ef_h = filter{1,1}; - ef_d = filter{1,2}; - ef_nh = filter{2,1}; - ef_nd = filter{2,2}; + ef_h = filters{1,1}; + ef_d = filters{1,2}; + ef_nh = filters{2,1}; + ef_nd = filters{2,2}; newSpikeTrain = [oldSpikeTrain(1:indx-1) timeToAdd oldSpikeTrain(indx:end)]; %possibly inefficient, change if problematic (only likely to be a problem for large numbers of spikes) diff --git a/utilities/get_initial_sample.m b/utilities/get_initial_sample.m index 285cbda..d1e0ae6 100644 --- a/utilities/get_initial_sample.m +++ b/utilities/get_initial_sample.m @@ -1,7 +1,16 @@ function SAM = get_initial_sample(Y,params) % obtain initial sample by performing sparse noise-constrained deconvolution -% Author: Eftychios A. Pnevmatikakis + +% Inputs: +% Y: observed fluorescence trace +% params: parameter struct (results of constrained foopsi can be passed here to avoid unnecessary computation) + +% Output: +% SAM: sample structure with values for spikes in continuous time, +% amplitude, baseline, noise, initial concentration and time constants + +% Author: Eftychios A. Pnevmatikakis, 2016, Simons Foundation if isfield(params,'p'); options.p = params.p; else options.p = 1; end diff --git a/utilities/get_next_spikes.m b/utilities/get_next_spikes.m index 9813975..39750df 100755 --- a/utilities/get_next_spikes.m +++ b/utilities/get_next_spikes.m @@ -1,193 +1,159 @@ -function [samples, ci] = get_next_spikes(curr_spikes,curr_calcium,calciumSignal,ef,tau,calciumNoiseVar, lam, proposalVar, add_move, Dt, A, con_lam) - - %addMoves, dropMoves, and timeMoves give acceptance probabilities for each subclass of move - %the samples will be a cell array of lists of spike times - the spike times won't be sorted but this shouldn't be a problem. - - %noise level here matters for the proposal distribution (how much it - %should trust data for proposals vs how much it should mix based on uniform prior) - %this is accounted for by calciumNoiseVar - +function [new_spikes, new_calcium, moves] = get_next_spikes(curr_spikes,curr_calcium,calciumSignal,filters,tau,calciumNoiseVar, lam, proposalSTD, add_move, Dt, A, con_lam) + +% Function for sampling the next of spike times given the current set of spike time and observed fluorescence trace + +% Inputs: +% curr_spikes: current set of spike times in continuous time +% curr_calcium: current estimate of noiseless calcium trace +% calciumSingal: observed fluorescence trace +% filters: exponential rise and decay kernels for calcium transient +% tau: continuous time rise and decay time constants +% calciumNoiseVar: observation noise variance +% lam: function handle for computing firing rate +% proposalSTD: standard deviation for perturbing a spike time +% add_move: # of addition/removal proposals per sample +% Dt: time-bin width +% A: spike amplitude +% con_lam: flag for constant firing rate (speeds up computation) + +% Outputs: +% new_spikes: new set of spike times +% new_calcium: new estimate of noiseless calcium trace +% moves: acceptance probabilities for perturbing,adding, droping spikes respectively + +% Author: Eftychios A. Pnevmatikakis and Josh Merel + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% initialize some parameters T = length(calciumSignal); %for all of this, units are bins and spiketrains go from 0 to T where T is number of bins ff = ~isnan(calciumSignal); % entries with data - -% nsweeps = 1e3; %number of sweeps. - nsweeps = 1; - samples = cell(nsweeps); - + %% start with initial spiketrain and initial predicted calcium si = curr_spikes; %initial set of spike times has no spikes - this will not be sorted but that shouldn't be a problem - ci = curr_calcium; %initial calcium is set to baseline + new_calcium = curr_calcium; %initial calcium is set to baseline N = length(si); %number of spikes in spiketrain %initial logC - compute likelihood initially completely - updates to likelihood will be local - logC = -norm(ci(ff)-calciumSignal(ff))^2; - - %m = p_spike*Dt*T; - + logC = -norm(new_calcium(ff)-calciumSignal(ff))^2; + %flag for uniform vs likelihood proposal (if using likelihood proposal, then time shifts are pure Gibbs) - %this really should be split into four cases, + %this really should be split into four cases (not implemented yet), % 1) RW for time shifts with uniform add/drop % 2) RW for time shifts with likelihood proposal add/drop % 3) Gibbs for time shifts with uniform add/drop % 4) Gibbs for time shifts with likelihood proposal add/drop %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - - - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - %% loop over sweeps to generate samples addMoves = [0 0]; %first elem is number successful, second is number total dropMoves = [0 0]; timeMoves = [0 0]; - for i = 1:nsweeps - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - %% loop over spikes, perform spike time move (could parallelize here for non-interacting spikes, i.e. spikes that are far enough away) + %% loop over spikes, perform spike time move (could parallelize here for non-interacting spikes, i.e. spikes that are far enough away) - %move - for ni = 1:N %possibly go through in a random order (if you like) - - tmpi = si(ni); - tmpi_ = si(ni)+(proposalVar*randn); %with bouncing off edges - if tmpi_<0 - tmpi_ = -(tmpi_); - elseif tmpi_>T - tmpi_ = T-(tmpi_-T); - end - - %set si_ to set of spikes with the move and ci_ to adjusted calcium and update logC_ to adjusted + for ni = 1:N %possibly go through in a random order (if you like) + + tmpi = si(ni); + tmpi_ = si(ni)+(proposalSTD*randn); %with bouncing off edges + if tmpi_<0 + tmpi_ = -(tmpi_); + elseif tmpi_>T + tmpi_ = T-(tmpi_-T); + end + + %set si_ to set of spikes with the move and ci_ to adjusted calcium and update logC_ to adjusted % [si_, ci_, logC_] = removeSpike(si,ci,logC,ef,tau,calciumSignal,tmpi,ni,Dt,A); % [si_, ci_, logC_] = addSpike(si_,ci_,logC_,ef,tau,calciumSignal,tmpi_,ni,Dt,A); - [si_, ci_, logC_] = replaceSpike(si,ci,logC,ef,tau,calciumSignal,tmpi,ni,tmpi_,Dt,A); - - %accept or reject - ratio = exp((logC_-logC)/(2*calciumNoiseVar)); - if ~con_lam; ratio = ratio*lam(tmpi)/lam(tmpi_); end - if ratio>1 %accept - si = si_; - ci = ci_; - logC = logC_; - timeMoves = timeMoves + [1 1]; - elseif rand1 %accept + si = si_; + new_calcium = ci_; + logC = logC_; + timeMoves = timeMoves + [1 1]; + elseif rand1 %accept + si = si_; + new_calcium = ci_; + logC = logC_; + addMoves = addMoves + [1 1]; + elseif rand0 + %propose a uniform removal + tmpi = randi(N); + [si_, ci_, logC_] = removeSpike(si,new_calcium,logC,filters,tau,calciumSignal,si(tmpi),tmpi,Dt,A); + + %reverse probability + rprob = 1/(T*Dt); + + %compute forward prob + fprob = 1/N; + %accept or reject - ratio = exp((1/(2*calciumNoiseVar))*(logC_ - logC))*(rprob/fprob)*lam(tmpi); %posterior times reverse prob/forward prob + ratio = exp((1/(2*calciumNoiseVar))*(logC_ - logC))*(rprob/fprob)*(1/lam(si(tmpi))); %posterior times reverse prob/forward prob + if ratio>1 %accept si = si_; - ci = ci_; + new_calcium = ci_; logC = logC_; - addMoves = addMoves + [1 1]; + dropMoves = dropMoves + [1 1]; elseif rand0 - %propose a uniform removal - tmpi = randi(N); - [si_, ci_, logC_] = removeSpike(si,ci,logC,ef,tau,calciumSignal,si(tmpi),tmpi,Dt,A); - - %reverse probability - rprob = 1/(T*Dt); - - %compute forward prob - fprob = 1/N; - - %accept or reject - ratio = exp((1/(2*calciumNoiseVar))*(logC_ - logC))*(rprob/fprob)*(1/lam(si(tmpi))); %posterior times reverse prob/forward prob - - if ratio>1 %accept - si = si_; - ci = ci_; - logC = logC_; - dropMoves = dropMoves + [1 1]; - elseif rand>>>>>> origin/master M = zeros(mS,nT); for i = 1:nT M(:,i) = hist(Mat(:,i),0:mS-1)/size(Mat,1); @@ -27,8 +38,8 @@ cmap = bone(100); cmap(2:26,:) = []; -fig = imagesc(M); axis xy; -<<<<<<< HEAD +fig = figure; +imagesc(M); axis xy; %set(gca,'Ytick',[0.5:(mS+.5)],'Yticklabel',[-1:(mS)]); %hold all; plot(traceData.spikeFrames+1); set(gca,'YLim',[1,5]) %set(gca,'YLim',[1.25,mS+.25]); hold all; scatter(1:nT,truespikes+1,[],'m'); @@ -37,27 +48,11 @@ set(gca,'Ytick',0.5+[-0.5:mS-.5],'Yticklabel',[-1+(0:mS)]); colormap(cmap); ylabel('# of Spikes ','fontweight','bold','fontsize',14); -%xlabel('Timestep ','fontweight','bold','fontsize',16); -title('Spike Histogram (MCMC) ','fontweight','bold','fontsize',14); -cbar = colorbar('Location','East'); -cpos = get(cbar,'Position'); -set(cbar,'Position',[(pos(1)+pos(3)-cpos(3)),cpos(2:4)]); -set(gca,'Xtick',[]) -set(cbar,'Color',[1,1,1]); -set(cbar,'Fontsize',12); -%xf = get(gca,'YTickLabel'); set(gca,'YTickLabel',xf,'fontsize',14) -======= -set(gca,'Ytick',[0.5:(mS-.5)],'Yticklabel',[-1:(mS-1)]); colormap(cmap); %hold all; plot(traceData.spikeFrames+1); set(gca,'YLim',[1,5]) -set(gca,'YLim',[1.5,mS]); -hold all; scatter(1:nT,truespikes+1,[],'m'); -set(gca,'YLim',[1.5,mS-1.5]); -pos = get(gca,'Position'); -set(gca,'Ytick',[-0.5:mS-1.5],'Yticklabel',[-1+(0:mS)]); -ylabel('# of Spikes ','fontweight','bold','fontsize',16); -%xlabel('Timestep ','fontweight','bold','fontsize',16); -title('Spike Histogram (MCMC) ','fontweight','bold','fontsize',16); +xlabel('Timestep ','fontweight','bold','fontsize',14); +title('Posterior Spike Histogram (MCMC) ','fontweight','bold','fontsize',14); cbar = colorbar('Location','East'); cpos = get(cbar,'Position'); -set(cbar,'Position',[(pos(1)+1.05*pos(3)),cpos(2:4)]); -set(gca,'Xtick',[]) ->>>>>>> origin/master +set(cbar,'Position',[pos(1)+pos(4),cpos(2:4)]); +%set(gca,'Xtick',[]) +%set(cbar,'Color',[1,1,1]); +set(cbar,'Fontsize',12); \ No newline at end of file diff --git a/utilities/removeSpike.m b/utilities/removeSpike.m index 3fb331c..bcc9f14 100755 --- a/utilities/removeSpike.m +++ b/utilities/removeSpike.m @@ -1,12 +1,33 @@ -function [newSpikeTrain, newCalcium, newLL] = removeSpike(oldSpikeTrain,oldCalcium,oldLL,filter,tau,obsCalcium,timeToRemove,indx,Dt,A) %#codegen - +function [newSpikeTrain, newCalcium, newLL] = removeSpike(oldSpikeTrain,oldCalcium,oldLL,filters,tau,obsCalcium,timeToRemove,indx,Dt,A) + +% Remove a given spike from the existing spike train. + +% Inputs: +% oldSpikeTrain: current spike train +% oldCalcium: current noiseless calcium trace +% oldLL: current value of the log-likelihood function +% filters: exponential rise and decay kernels for calcium transient +% tau: continuous time rise and decay time constants +% obsCalcium: observed fluorescence trace +% timetoRemove: time of the spike to be removed +% indx: place where the spike to be removed is in the existing spike train vector +% Dt: time-bin width +% A: spike amplitude + +% Outputs: +% newSpikeTrain: new vector of spike times +% newCalcium: new noiseless calcium trace +% newLL: new value of the log-likelihood function + +% Author: Eftychios A. Pnevmatikakis and Josh Merel + tau_h = tau(1); tau_d = tau(2); - ef_h = filter{1,1}; - ef_d = filter{1,2}; - ef_nh = filter{2,1}; - ef_nd = filter{2,2}; + ef_h = filters{1,1}; + ef_d = filters{1,2}; + ef_nh = filters{2,1}; + ef_nd = filters{2,2}; newSpikeTrain = oldSpikeTrain; newSpikeTrain(indx) = []; diff --git a/utilities/replaceSpike.m b/utilities/replaceSpike.m index f962c66..8d17676 100644 --- a/utilities/replaceSpike.m +++ b/utilities/replaceSpike.m @@ -1,5 +1,27 @@ function [newSpikeTrain, newCalcium, newLL] = replaceSpike(oldSpikeTrain,oldCalcium,oldLL,filter,tau,obsCalcium,timetoRemove,indx,timetoAdd,Dt,A) %#codegen - + +% Replace a given spike with a new one in the existing spike train. + +% Inputs: +% oldSpikeTrain: current spike train +% oldCalcium: current noiseless calcium trace +% oldLL: current value of the log-likelihood function +% filters: exponential rise and decay kernels for calcium transient +% tau: continuous time rise and decay time constants +% obsCalcium: observed fluorescence trace +% timetoRemove: time of the spike to be removed +% indx: place where the spike to be removed is in the existing spike train vector +% timetoAdd: time of the spike to be added +% Dt: time-bin width +% A: spike amplitude + +% Outputs: +% newSpikeTrain: new vector of spike times +% newCalcium: new noiseless calcium trace +% newLL: new value of the log-likelihood function + +% Author: Eftychios A. Pnevmatikakis and Josh Merel + tau_h = tau(1); tau_d = tau(2); @@ -33,7 +55,6 @@ end newCalcium(tmp) = newCalcium(tmp) - wef_h; end - %%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%% %handle ef_d next diff --git a/utilities/sampleSpikes.m b/utilities/sampleSpikes.m deleted file mode 100755 index d557dc5..0000000 --- a/utilities/sampleSpikes.m +++ /dev/null @@ -1,192 +0,0 @@ -function [samples, addMoves, dropMoves, timeMoves, N_sto] = sampleSpikes(calciumSignal,ef,tau,b,calciumNoiseVar, p_spike, proposalVar, nsweeps) - - %addMoves, dropMoves, and timeMoves give acceptance probabilities for each subclass of move - %the samples will be a cell array of lists of spike times - the spike times won't be sorted but this shouldn't be a problem. - - %noise level here matters for the proposal distribution (how much it - %should trust data for proposals vs how much it should mix based on uniform prior) - %this is accounted for by calciumNoiseVar - - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - %% initialize some parameters - T = length(calciumSignal); %for all of this, units are bins and spiketrains go from 0 to T where T is number of bins - -% nsweeps = 1e3; %number of sweeps. - samples = cell(nsweeps); - N_sto = []; - objective = []; - - %% start with initial spiketrain and initial predicted calcium - si = []; %initial set of spike times has no spikes - this will not be sorted but that shouldn't be a problem - ci = b*ones(1,T); %initial calcium is set to baseline - - N = length(si); %number of spikes in spiketrain - - %initial logC - compute likelihood initially completely - updates to likelihood will be local - logC = -(ci-calciumSignal)*(ci-calciumSignal)'; - - m = p_spike*T; - - %flag for uniform vs likelihood proposal (if using likelihood proposal, then time shifts are pure Gibbs) - %this really should be split into four cases, - % 1) RW for time shifts with uniform add/drop - % 2) RW for time shifts with likelihood proposal add/drop - % 3) Gibbs for time shifts with uniform add/drop - % 4) Gibbs for time shifts with likelihood proposal add/drop - - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - - - - - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - %% loop over sweeps to generate samples - addMoves = [0 0]; %first elem is number successful, second is number total - dropMoves = [0 0]; - timeMoves = [0 0]; - for i = 1:nsweeps - - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - %% loop over spikes, perform spike time move (could parallelize here for non-interacting spikes, i.e. spikes that are far enough away) - - %move - for ni = 1:N %possibly go through in a random order (if you like) - - tmpi = si(ni); - tmpi_ = si(ni)+(proposalVar*randn); %with bouncing off edges - if tmpi_<0 - tmpi_ = -(tmpi_); - elseif tmpi_>T - tmpi_ = T-(tmpi_-T); - end - - %set si_ to set of spikes with the move and ci_ to adjusted calcium and update logC_ to adjusted - [si_, ci_, logC_] = removeSpike(si,ci,logC,ef,tau,calciumSignal,tmpi,ni); - [si_, ci_, logC_] = addSpike(si_,ci_,logC_,ef,tau,calciumSignal,tmpi_); - - %accept or reject - ratio = exp((1/(2*calciumNoiseVar))*(logC_-logC)); - if ratio>1 %accept - si = si_; - ci = ci_; - logC = logC_; - timeMoves = timeMoves + [1 1]; - elseif rand1 %accept - si = si_; - ci = ci_; - logC = logC_; - addMoves = addMoves + [1 1]; - elseif rand0 - %propose a uniform removal - tmpi = randi(N); - [si_ ci_, logC_] = removeSpike(si,ci,logC,ef,tau,calciumSignal,si(tmpi),tmpi); - - %reverse probability - rprob = 1/T; - - %compute forward prob - fprob = 1/N; - - %accept or reject - ratio = exp((1/(2*calciumNoiseVar))*(logC_ - logC))*(rprob/fprob)*((T-m)/m); %posterior times reverse prob/forward prob - - if ratio>1 %accept - si = si_; - ci = ci_; - logC = logC_; - dropMoves = dropMoves + [1 1]; - elseif rand