how to avoid NaN in optim() - R-help

Fri, Oct 1, 2010 11:22 AM #

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

 ____________________________________________________________________
 
 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.