"By popular demand": likelihood [🇷 for BE/BA]

posted by ElMaestro  – Denmark, 2020-07-24 10:07 (792 d 15:16 ago) – Posting: # 21783
Views: 12,959

Hi all,

to make the objective function return the likelihood, use e.g.:

Obj.F12=function(Pars)
##returns the restricted log likelihood.
##apart from k, this is line 3 of page 10 of:
##http://people.csail.mit.edu/xiuming/docs/tutorials/reml.pdf
##k is invariant with parameter estimates and variance components.
{
  CovM=Create.CovM(Pars)
  k=  -((length(y)-ncol(X))/2)*log(2*pi)
  A= -0.5*log( det(CovM))
  B= -0.5*log( det(t(X) %*% solve(CovM) %*% X))
  est.b = solve(t(X) %*% solve(CovM) %*% X) %*% t(X) %*% solve(CovM) %*% y
  tmp= y - X %*% est.b
  C=-0.5 *(t(tmp) %*% solve(CovM) %*% tmp)
  return(k+A+B+C)
}


The optimised object (F in my original code) will hold the optimised likelihood as F$value.
Opportunity: Will run a tiny bit faster if k is public, and if solve(CovM) is done a single time.


For example at reltol=1e-8 I get with EMA dataset II a value of 15.33728.

This matches exactly what SAS reports, see post # 21670 above.

Pass or fail!
ElMaestro

Complete thread:

UA Flag
Activity
 Admin contact
22,385 posts in 4,684 threads, 1,594 registered users;
online 11 (0 registered, 11 guests [including 10 identified bots]).
Forum time: Sunday 01:24 CEST (Europe/Vienna)

You really don’t know what you don’t know until you write about it.
Then, everyone knows what you don’t know.    Rod Machado

The Bioequivalence and Bioavailability Forum is hosted by
BEBAC Ing. Helmut Schütz
HTML5