New script:
load("Temperature.mat")
T = T'; % from column vector to row vector
AT = T-min(T); % this is the adjusted positive temperatures
r = mypoiss(5,0.5,1,AT);
s = 0:3; % lag range
om = 0.1:0.1:0.4; % decay rate range
mu = 0.1:0.1:0.4; % baseline rate range
Ns = numel(s);
Nom = numel(om);
Nmu = numel(mu);
nll = nan(Ns,Nom,Nmu);
for si = 1:Ns
omi = 1:Nom;
mui = 1:Nmu;
nll(si,omi,mui) = mylikelihood2(s(si),om,mu,AT,r,Nom,Nmu);
end
[ss,omm,muu]=meshgrid(s,om,mu);
scatter3(ss(:),omm(:),muu(:),[],nll(:))
New function:
function nll = mylikelihood2(s,om,mu,AT,r,Nom,Nmu)
% output; nll is the negative log likelihood
% input; s is lag, om is decay rate, mu is baseline rate
% T is temperatures, t = time in days
t = max(s)+1:numel(AT);
[omm,muu] = meshgrid(om,mu);
nll = zeros(Nom,Nmu);
for ti = 1:numel(t)
k = abs(s:-1:(s-t(ti)+1));
id = ((t(ti)-1):-1:0)+1;
lambda = muu+sum(omm(:).^k.*AT(id))/sum(omm(:).^k);
% store all the results
nll = nll - (r(ti).*log(lambda)) + (lambda);
end
