[R] how to avoid NaN in optim()

[Ravi Varadhan]

Change the `else NA' to `else Inf' in your loglik function  and then run the following:


library(BB)

p0 <- runif(2)

spg(p0, lik (176,182 , 60 ,17) , lower=0,  upper=1)


This will give you correct results.

Ravi.

____________________________________________________________________

Ravi Varadhan, Ph.D.
Assistant Professor,
Division of Geriatric Medicine and Gerontology
School of Medicine
Johns Hopkins University

Ph. (410) 502-2619
email: rvaradhan at jhmi.edu


----- Original Message -----
From: arindam fadikar <arindam.fadikar at gmail.com>
Date: Thursday, September 30, 2010 2:17 pm
Subject: [R] how to avoid NaN in optim()
To: r-help at r-project.org


> hi ,
>  
>  lik <- function(nO, nA, nB, nAB){
>  
>  loglik <- function(par)
>    {
>        p=par[1]
>        q=par[2]
>        r <- 1 - p - q
>  
>        if (c(p,q,r) > rep(0,3) && c(p,q,r) < rep(1,3) )
>  
>          {
>            -(2 * nO * log (r) + nA * log (p^2 + 2 * p * r)
>              + nB * log (q^2 + 2 * q * r)
>              + nAB * (log(2) +log(p) +log(q)))
>         }
>        else
>          NA
>    }
>  
>  loglik
>  
>  }
>  
>  
>  i want to maximize this likelihood function over the range (0,0,0) to
>  (1,1,1). so i tried
>  
>  optim ( c(0.3,0.3), lik ( 176,182 , 60 ,17) , method = "CG")
>  
>  and obtained the following :
>  
>  > optim(c(0.3,0.3), fn, method="CG")
>  $par
>  [1] 0.26444187 0.09316946
>  
>  $value
>  [1] 492.5353
>  
>  $counts
>  function gradient
>       130       19
>  
>  $convergence
>  [1] 0
>  
>  $message
>  NULL
>  
>  Warning messages:
>  1: In log(q^2 + 2 * q * r) : NaNs produced
>  2: In log(q) : NaNs produced
>  3: In log(q^2 + 2 * q * r) : NaNs produced
>  4: In log(q) : NaNs produced
>  
>  
>  please help ...
>  
>  
>  -- 
>  Arindam Fadikar
>  M.Stat
>  Indian Statistical Institute.
>  New Delhi, India
>  
>  	[[alternative HTML version deleted]]
>  
>  ______________________________________________
>  R-help at r-project.org mailing list
>  
>  PLEASE do read the posting guide 
>  and provide commented, minimal, self-contained, reproducible code.