1;

## Calculate the log of the P(Y) probability for one subject. 
function v = loglike1(Y, sig2, mu, tau2, T)

  ## The lambda grid points.
  tau = sqrt(tau2);
  h = 6*tau / 20;
  L = [mu-3*tau:h:mu+3*tau]';

  ## Cycle through the lambda grid points.
  ep = zeros(length(L),1);
  for k=1:length(L)

    ## The expected values of the Y components.
    E = exp(L(k)*T);

    ## The residuals.
    R = Y - E;

    ## The exponent.
    ep(k) = -0.5*sumsq(R)/sig2 - 0.5*(L(k)-mu)^2/tau2;

  endfor

  epmax = max(ep);
  ep = ep - epmax;

  ## Quadrature weights.
  W = ones(length(L),1);
  W(1) = 1/2;
  W(length(W)) = 1/2;
  W = W*h;

  v = epmax + log(sum(W.*exp(ep))) - log(tau2)/2 - size(Y,2)*log(sig2)/2; 

endfunction


## Calculate the log of the P(Y) probability for all subjects. 
function v = loglike(Y, sig2, mu, tau2, T)

  v = 0;
  
  for i=1:size(Y,1)
    v = v + loglike1(Y(i,:)', sig2, mu, tau2, T);
  endfor

endfunction


## Calculate the MLE of mu using bisection.
function mu = est_mu(Y, sig2, tau2, T)

  n = size(Y,1);
  q = size(Y,2);

  ## A crude initial estimate for the center of the bracket.
  YL = log(Y(:,q)) / T(q);
  YL = YL(find(Y(:,q) > 1e-4));
  mu2 = mean(real(YL));
  f2 = loglike(Y, sig2, mu2, tau2, T);

  ## Get a left bracket point.
  mu1 = mu2-1;
  while (1)
    f1 = loglike(Y, sig2, mu1, tau2, T);
    if (f1 < f2)
      break;
    endif
    mu1 = mu1-1;
  endwhile

  ## Get a right bracket point.
  mu3 = mu2+1;
  while (1)
    f3 = loglike(Y, sig2, mu3, tau2, T);
    if (f3 < f2)
      break;
    endif
    mu3 = mu3+1;
  endwhile

  ## Squeeze the bracket.
  while (mu3-mu1 > 1e-4)

    if (mu3-mu2 > mu2-mu1)

      mm = (mu2+mu3)/2;
      ff = loglike(Y, sig2, mm, tau2, T);
      if (ff > f2)
	mu1 = mu2;
	f1 = f2;
	mu2 = mm;
	f2 = ff;
      else
	mu3 = mm;
	f3 = ff;
      endif

    else

      mm = (mu1+mu2)/2;
      ff = loglike(Y, sig2, mm, tau2, T);
      if (ff > f2)
	mu3 = mu2;
	f3 = f2;
	mu2 = mm;
	f2 = ff;
      else
	mu1 = mm;
	f1 = ff;
      endif

    endif

  endwhile

  mu = mu2;

endfunction



## Simulate data from the population.
function Y = simdat(sig2, mu, tau2, n, T)

  L = sqrt(tau2)*randn(n,1) + mu;
  EY = e.^(L * T');
  Y = EY + sqrt(sig2)*randn(size(EY));

endfunction

## Population structure.
mu = -1;
sig2 = 0.3;
tau2 = 0.5;
T = [0.1:0.1:1]';
n = 300;

## Run the simulation study.
for r=1:100
  Y = simdat(sig2, mu, tau2, n, T);
  mu_mle(r) = est_mu(Y, sig2, tau2, T);
  disp(mu_mle(r));
endfor