% This demonstrates a non-linear regression model for the onset of rigor mortis. 

% Niderkorn's (1872) observations on 113 bodies provides the main reference
% database for the development of rigor mortis and is commonly
% cited in textbooks. His data was as follows (Ref. 19 at p. 31): 
	      
% Here's the data:
bodies= [2 14 31 14 20 11 7 4 7 1 1 2]
hours=[2 3 4 5 6 7 8 9 10 11 12 13]
% we compute the cumulative number of deaths by hour:
n=size(hours)(2);
cs(1)=bodies(1);
for i=2:n
  cs(i)=cs(i-1)+bodies(i);
endfor
cs
    
% Having examined the data, here's the model that we believe fits the data:
function f = f(x, gamma, alpha, beta)
  f=gamma * exp(- (alpha / (x ^ beta)));
endfunction

% one can transform the data and use linear regression to get a good start for the model
% (if we can use a guess of our own for gamma -- here 120):
gamma=120 % this was our guess for the horizontal asymptote, based on the data

A=[ones(n,1) ln(hours)'];
b=ln(- ln(cs/120))';
reg=A\b
xstar=A*reg;

% Then we deduce the other two parameters:
alpha=exp(reg(1))
beta=-reg(2)

% and we compute the model using the results obtained from the transformed
% linear regression:
t = 0.01:.05:14;
ft=gamma * exp(- (alpha ./ (t .^ beta)));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
% Here's the non-linear regression routine:
% 
% First, we'll need some partial derivatives:
function fpg = fpg(x, gamma, alpha, beta)
  fpg=exp (- (alpha / (x ^ beta)));
endfunction
function fpa = fpa(x, gamma, alpha, beta)
  fpa=gamma * exp (- (alpha / (x ^ beta))) * (- (1 / (x ^ beta)));
endfunction
function fpb = fpb(x, gamma, alpha, beta)
  fpb=gamma * exp (- (alpha / (x ^ beta))) * (- ((- (alpha * ((exp (beta * log(x))) * log(x)))) / ((x ^ beta) ^ 2)));
endfunction

% Now the iterative scheme that we'll do again and again, until satisfied:
function nlr = nlr(runs, hours, cs, bodies,gamma,alpha,beta)
  n=size(hours)(2);
  for run= 1:runs
	      % Let's compute our estimates, and the partial derivatives:
    for hour=1:n
      estimate(hour)=f(hours(hour),gamma,alpha,beta);
      fpgs(hour)=fpg(hours(hour),gamma,alpha,beta);
      fpas(hour)=fpa(hours(hour),gamma,alpha,beta);
      fpbs(hour)=fpb(hours(hour),gamma,alpha,beta);
    endfor
	      % perform another linear regression:
    A=[fpgs' fpas' fpbs'];
    b=cs' - estimate';
    reg=A\b;
	      % update the coefficients:
    coefs= [gamma alpha beta] + reg';
    gamma=coefs(1)
    alpha=coefs(2)
    beta=coefs(3)
  endfor
  nlr=[gamma alpha beta];
endfunction
    
% Run the non-linear regression four times and get the updated parameter estimates:
% We use the linear transformed values for starting values
coefs=nlr(4, hours, cs, bodies,gamma,alpha,beta);
gamma=coefs(1)
alpha=coefs(2)
beta=coefs(3)

% now generate a vector of model predictions:
ftupdate=gamma * exp(- (alpha ./ (t .^ beta)));

% Let's plot the results: first of all, the linear regression picture:
subplot(2,2,1)
xlabel("hours since death")
ylabel("cumulative bodies in rigor mortis")
title("Niderkorn's (1872), 113 bodies")
plot(hours,cs,"3*")
clg

% Let's plot the results: first of all, the linear regression picture:
subplot(2,2,2)
xlabel("ln(hours) since death")
ylabel("ln(- ln(cumulative/120)) bodies")
title("transformed data for the onset of rigor mortis")
plot(ln(hours),b,"3*")
clg

% Let's plot the results: first of all, the linear regression picture:
subplot(2,2,3)
xlabel("ln(hours) since death")
ylabel("ln(- ln(cumulative)) bodies")
title("Linear regression model (transformed)")
plot(ln(hours),b,"3*;b;",ln(hours),xstar,"x;x*;",ln(hours),xstar,"3-;model;")
clg

% And now the results of the non-linear regression(s):
subplot(2,2,4)
xlabel("hours since death")
ylabel("proportion in rigor mortis")
title("Non-linear regression model")
plot(hours,cs/gamma,"3*;cumulative;",t,ft/gamma,"3-;transformed linear model;",t,ftupdate/gamma,"2-;non-linear model;")