0001 function delay = PETHTransition(psth,timeBins,varargin)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036
0037
0038
0039
0040
0041
0042
0043
0044 method = 'ml';
0045 show = false;
0046
0047
0048 if nargin < 2 | mod(length(varargin),2) ~= 0,
0049 error('Incorrect number of parameters (type ''help <a href="matlab:help PETHTransition">PETHTransition</a>'' for details).');
0050 end
0051
0052
0053 for i = 1:2:length(varargin),
0054 if ~ischar(varargin{i}),
0055 error(['Parameter ' num2str(i+1) ' is not a property (type ''help <a href="matlab:help PETHTransition">PETHTransition</a>'' for details).']);
0056 end
0057 switch(lower(varargin{i})),
0058
0059 case 'method',
0060 method = lower(varargin{i+1});
0061 if ~isastring(method,'ml','ls','pl'),
0062 error('Incorrect value for property ''method'' (type ''help <a href="matlab:help PETHTransition">PETHTransition</a>'' for details).');
0063 end
0064
0065 case {'show','plot'},
0066 show = lower(varargin{i+1});
0067 if ~isastring(show,'on','off'),
0068 error('Incorrect value for property ''show'' (type ''help <a href="matlab:help PETHTransition">PETHTransition</a>'' for details).');
0069 end
0070
0071 otherwise,
0072 error(['Unknown property ''' num2str(varargin{i}) ''' (type ''help <a href="matlab:help PETHTransition">PETHTransition</a>'' for details).']);
0073
0074 end
0075 end
0076
0077 if strcmp(method,'ml'),
0078
0079
0080
0081 f = psth;
0082 n = length(f);
0083 for i = 1:n,
0084 F(i) = sum(log(1:max([1 f(i)])));
0085 end
0086 for theta = 1:n-1,
0087 t1 = 1:theta;t1 = t1';
0088 t2 = theta+1:n;t2 = t2';
0089 lambda1 = sum(f(t1))/(theta+1);
0090 lambda2 = sum(f(t2))/(n-theta);
0091 likelihood(theta) = ...
0092 - lambda1 * (theta+1) ...
0093 + log(lambda1) * sum(f(t1)) ...
0094 - sum(F(t1)) ...
0095 - lambda2 * (n-theta) ...
0096 + log(lambda2) * sum(f(t2)) ...
0097 - sum(F(t2));
0098 end
0099 likelihood(end+1) = likelihood(end);
0100 d = find(likelihood == max(likelihood));
0101 d = d(1);
0102 delay = timeBins(d);
0103
0104 if show,
0105 fig = figure;
0106 set(fig,'name','PETH Transition - Maximum Likelihood','number','off');
0107 subplot(2,1,1);hold on;
0108 bar(timeBins,psth);
0109 plot([delay delay],ylim,'k','linestyle','--');
0110 title(['delay = ' num2str(delay) ' s']);
0111 ylabel('Occurrences');
0112 subplot(2,1,2);hold on;
0113 plot(timeBins,likelihood,'Color',[1 0 0],'Marker','none','LineStyle','-');
0114 plot([delay delay],ylim,'k','linestyle','--');
0115 xlabel('Time (s)');
0116 ylabel('Likelihood');
0117 end
0118
0119 elseif strcmp(method,'ls'),
0120
0121
0122
0123 f0 = psth;
0124 n0 = length(f0);
0125 t0 = 1:n0;t0 = t0';
0126 for i = 1:n0,
0127 F0(i) = sum(f0(1:i-1));
0128 end
0129 F0 = F0';
0130 start = 1;
0131 stop = length(f0);
0132 f = f0(start:stop);
0133 n = length(f);
0134 for i = 1:length(f),
0135 F(i) = sum(f(1:i-1));
0136 end
0137 F = F';
0138 t = 1:n;t = t';
0139 squaredSum(1:n) = NaN;
0140 for theta = 2:n-1,
0141 t1 = 1:theta; t1 = t1';
0142 t2 = theta+1:n; t2 = t2';
0143 lambda1(theta) = ...
0144 (...
0145 - sum(t1 .* F(t1)) ...
0146 - theta * sum(F(t2)) ...
0147 + theta * sum(theta-t2) * sum(F(t2) .* (theta-t2)) / sum((theta-t2).^2) ...
0148 ) / (...
0149 + theta^2 * sum(theta-t2) * sum(theta-t2) / sum((theta-t2).^2) ...
0150 - sum(t1.^2) ...
0151 - theta^2 * (n-theta) ...
0152 );
0153 lambda2(theta) = ...
0154 (...
0155 + lambda1(theta) * theta * sum(theta-t2) ...
0156 - sum(F(t2).*(theta-t2)) ...
0157 ) ...
0158 / sum((theta-t2).^2);
0159 squaredSum(theta) = ...
0160 + sum((F(t1)-lambda1(theta)*t1).^2) ...
0161 + sum((F(t2)-lambda1(theta)*theta+lambda2(theta)*(theta-t2)).^2);
0162 end
0163 d = find(squaredSum == min(squaredSum));
0164 d = d(1);
0165 delay = timeBins(d);
0166
0167 if show,
0168 fig = figure;
0169 set(fig,'name','PETH Transition - Least Square Error','number','off');
0170 subplot(2,1,1);hold on;
0171 bar(timeBins,psth);
0172 plot([delay delay],ylim,'k','linestyle','--');
0173 title(['delay = ' num2str(delay) ' s']);
0174 ylabel('Occurrences');
0175 subplot(2,1,2);hold on;
0176 t1 = 1:d;t1 = t1';
0177 t2 = d+1:n;t2 = t2';
0178 plot(timeBins(t0),F0,'Color',[1 0 0],'Marker','none','LineStyle','-');
0179 plot(timeBins(t1+start-1),lambda1(d)*(t1)+F0(start),'Color',[0 0 0]);
0180 plot(timeBins(t2+start-1),lambda2(d)*(t2-d)+lambda1(d)*(d)+F0(start),'Color',[0 0 0]);
0181 plot([delay delay],ylim,'k','linestyle','--');
0182 xlabel('Time (s)');
0183 ylabel('Cumulative count');
0184 end
0185
0186 elseif strcmp(method,'pl')
0187
0188
0189
0190 f = psth;
0191
0192 n = length(f);
0193 t = 1:n;t = t';
0194 for theta = 2:n-2,
0195 t1 = 1:theta;t1 = t1';
0196 t2 = theta+1:n;t2 = t2';
0197 a1(theta) = LinearRegression(t1,f(t1));
0198 a2(theta) = LinearRegression(t2,f(t2),a1(theta)*theta);
0199 squaredSum(theta) = sum((f(t1)-a1(theta)*t1).^2) ...
0200 + sum((f(t2)-a2(theta)*t2-a1(theta)*theta).^2);
0201 end
0202 squaredSum(1) = Inf;
0203 d = find(squaredSum == min(squaredSum));
0204 d = d(1);
0205 delay = timeBins(d);
0206
0207 if show,
0208 fig = figure;
0209 set(fig,'name','PETH Transition - Piecewise Linear Fit','number','off');
0210 subplot(2,1,1);hold on;
0211 bar(timeBins,psth);
0212 plot([delay delay],ylim,'k','linestyle','--');
0213 title(['delay = ' num2str(delay) ' s']);
0214 ylabel('Occurrences');
0215 subplot(2,1,2);hold on;
0216 t1 = 1:d;t1 = t1';
0217 t2 = d+1:n;t2 = t2';
0218 plot(timeBins(t),f,'Color',[1 0 0],'Marker','o','LineStyle','none');
0219 plot(timeBins(t1),a1(d)*t1,'k');
0220 plot(timeBins(t2),a2(d)*t2+a1(d)*d,'k');
0221 plot([delay delay],ylim,'k','linestyle','--');
0222 xlabel('Time (s)');
0223 ylabel('Occurrences');
0224 end
0225
0226 end
0227
0228 function [a,b] = LinearRegression(x,y,b)
0229
0230 if nargin == 2, b = 0; end
0231
0232 y = y - b;
0233 X = sum(x);
0234 X2 = sum(x.^2);
0235 Y = sum(y);
0236 XY = sum(x.*y);
0237 n = length(x);
0238
0239 b = (Y*X2-X*XY) / (n*X2-X.^2);
0240 a = (n*XY-X*Y) / (n*X2-X.^2);