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

posted by ElMaestro  – Denmark, 2020-07-24 10:07 (672 d 11:39 ago) – Posting: # 21783
Views: 12,553

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,110 posts in 4,630 threads, 1,567 registered users;
online 15 (0 registered, 15 guests [including 2 identified bots]).
Forum time: Friday 21:47 CEST (Europe/Vienna)

We absolutely must leave room for doubt
or there is no progress and no learning.
There is no learning without having to pose a question.
And a question requires doubt.
People search for certainty.
But there is no certainty.    Richard Feynman

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