https://github.com/hasselmonians/mle_rhythmicity
Tip revision: 7c50ba73fa34dbe2586fb0f9e0204bd1e28f3209 authored by jrclimer on 11 July 2016, 16:29:52 UTC
Resolves #7
Resolves #7
Tip revision: 7c50ba7
cif_generator.m
function [ cif, cif_int ] = cif_generator( funname )
%CIF_GENERATOR Generate function handles for rhythm conditional intensity
% This helper function generates function handles for fitting rhythmic
% distributions of lags. Using functional programming in this way greatly
% speeds up evaluation during runtime.
%
% INPUT
% funname: String name of the CIF handles to be generated. Can be:
% flat - Exponential decay with no oscillation
% pure - Pure sinusoidal decay
% noskip - Decaying rhythm with no skipping
% full - Decaying rhythm with skipping
%
% OUTPUT
% cif - Function handle for the Conditional Intensity Function (CIF).
% with the form @(lags_list,...). lags_list is the lags at which
% to evaluate the CIF, and it is followed by the appropriate
% parameters
% cif_int - The function handle for the definite integral of the CIF,
% with the form @(max_lag,...).
%
% Copyright 2015-2016 Trustees of Boston University
% All rights reserved.
%
% This file is part of mle_rhythmicity revision 2.0. The last committed
% version of the previous revision is the SHA starting with 93862ac...
%
% This code has been freely distributed by the authors under the BSD
% license (http://opensource.org/licenses/BSD2-Clause). If used or
% modified, we would appreciate if you cited our papers:
%
% Climer JR, DiTullio R, Newman EL, Hasselmo ME, Eden UT. (2014),
% Examination of rhythmicity of extracellularly recorded neurons in the
% entorhinal cortex. Hippocampus, 25:460-473. doi: 10.1002/hipo.22383.
%
% Hinman et al., Multiple Running Speed Signals in Medial Entorhinal
% Cortex, Neuron (2016). http://dx.doi.org/10.1016/j.neuron.2016.06.027
switch funname
case 'flat'% Non-rhythmic
cif = @(lags_list,tau,b)(1-b)*exp(-lags_list*10^-tau)+b;
cif_int = @(max_lag,tau,b)10^tau*(1-b)*(1-exp(-max_lag*10^-tau))+b*max_lag;
case 'pure'% Pure sinusoid
cif = @(lags_list,f)cos(2*pi*f*lags_list)+1;
cif_int = @(max_lag,f)sin(2*pi*f*max_lag)/(2*pi*f)+max_lag;
case 'noskip'% Non-skipping
cif = @(lags_list,tau,b,c,f,r)(1-b).*exp(-lags_list.*10.^-tau).*(r.*exp(-lags_list.*10.^-c).*...
cos(2*pi*f.*lags_list)...
+1)+b;
cif_int = @(max_lag,tau,b,c,f,r)b*max_lag+(1-b).*10.^tau.*(1-exp(-max_lag*10.^-tau))+r.*(1-b).*...
exp(-max_lag*(10.^-c+10.^-tau)).*((10.^-c+10.^-tau).*(exp(max_lag.*(10.^-c+10.^-tau))-cos(2*pi*f*max_lag))+2*pi*f.*sin(2*pi*f*max_lag))./...
(4*pi^2*f.^2+(10.^-c+10.^-tau).^2);
case 'full'
fixb = @(b)max(b,realmin);% b>0
cif = @(lags_list,tau,b,c,f,s,r)(1-b).*exp(-lags_list.*10.^-tau).*(r.*exp(-lags_list.*10.^-c).*...
((2+2*sqrt(1-s)-s).*cos(2*pi*f.*lags_list)+4*s.*cos(pi*f.*lags_list)+2-2*sqrt(1-s)-3*s)/4 ...
+1)+b;
cif = @(lags_list,tau,b,c,f,s,r)cif(lags_list,tau,fixb(b),c,f,s,r);
cif_int = @skipping_integral;
otherwise
end
end
% Integral of the full distribution with skipping
function D = skipping_integral(max_lag,tau,b,c,f,s,r)
A = 10.^-c+10.^-tau;
B = exp(max_lag*A);
b = max(b,realmin);
D=b*max_lag+...
...
(1-b).*...
(10.^tau).*...
(1-exp(-max_lag*10.^-tau))+...
(r.*(1-b)/4).*...
((2+2*sqrt(1-s)-s).*((B.^-1).*(A.*(B-cos(2*pi*max_lag*f))+2*pi*f.*sin(2*pi*max_lag*f))./...
(4*pi^2*f.^2+A.^2))+...
4*s.*((B.^-1).*(A.*(B-cos(pi*max_lag*f))+pi*f.*sin(pi*max_lag*f))./...
(pi^2*f.^2+A.^2))+...
(2-2*sqrt(1-s)-3*s).*(10.^(c+tau)).*(1-(B.^-1))./(10.^c+10.^tau));
end