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

posted by ElMaestro  – Belgium?, 2020-07-24 10:07 (68 d 21:31 ago) – Posting: # 21783
Views: 7,461

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.

I could be wrong, but...

Best regards,
ElMaestro

R's base package has 274 reserved words and operators, along with 1761 functions. I can use 18 of them (about 14 of them properly). I believe this makes me the Donald Trump of programming.

Complete thread:

Activity
 Admin contact
21,090 posts in 4,398 threads, 1,469 registered users;
online 18 (0 registered, 18 guests [including 9 identified bots]).
Forum time: Thursday 07:38 CEST (Europe/Vienna)

In these days, a man who says a thing cannot be done
is quite apt to be interrupted by some idiot doing it.    Elbert Green Hubbard

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