[R] how to avoid NaN in optim()

[Ravi Varadhan]

You do not have to worry about re-defining your function to handle the
constraints on parameters.  You can let the optimizer worry about it.

lik <- function(nO, nA, nB, nAB){
	loglik <- function(par)
	{
	p=par[1]
	q=par[2]
	r <- 1 - p - q

	-(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)))
	}
}


library(BB)
# constraints:  x[1] > 0; x[2] > 0; x[1] + x[2] < 1

Amat <- matrix(c(1,0,0,1,-1,-1), 3, 2, byrow=TRUE)

b <- c(0, 0, -1)

c1 <- runif(1) 

p0 <- c(c1, max(0.01, 0.99-c1)) 

spg(p0, lik ( 176,182,60,17) , project="projectLinear",
projectArgs=list(A=Amat, b=b, meq=0))

Hope this helps,
Ravi.

-----Original Message-----
From: Ravi Varadhan [mailto:rvaradhan at jhmi.edu] 
Sent: Thursday, September 30, 2010 3:09 PM
To: Ravi Varadhan
Cc: arindam fadikar; r-help at r-project.org
Subject: Re: [R] how to avoid NaN in optim()

Here is how you do it:

library(BB)

Amat <- matrix(c(1,0,0,1,-1,-1), 3, 2, byrow=TRUE)

b <- c(0, 0, -1)

p0 <- c(0.5, 0.4)

spg(p0, lik ( 176,182 , 60 ,17) , project="projectLinear",
projectArgs=list(A=Amat, b=b, meq=0))


Hope this helps,
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: Ravi Varadhan <rvaradhan at jhmi.edu>
Date: Thursday, September 30, 2010 3:04 pm
Subject: Re: [R] how to avoid NaN in optim()
To: Ravi Varadhan <rvaradhan at jhmi.edu>
Cc: arindam fadikar <arindam.fadikar at gmail.com>, r-help at r-project.org


> You also need the constrain that par[1] + par[2] < 1 in order to avoid 
> NaNs.
>  
>  You can do this using the `projectLinear' argument in  `spg'.
>  
>  library(BB)
>  ?spg
>  
>  
>  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: Ravi Varadhan <rvaradhan at jhmi.edu>
>  Date: Thursday, September 30, 2010 2:54 pm
>  Subject: Re: [R] how to avoid NaN in optim()
>  To: arindam fadikar <arindam.fadikar at gmail.com>
>  Cc: r-help at r-project.org
>  
>  
>  > 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.
>  >  
>  >  ______________________________________________
>  >  R-help at r-project.org mailing list
>  >  
>  >  PLEASE do read the posting guide 
>  >  and provide commented, minimal, self-contained, reproducible code. 
>