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:rvarad...@jhmi.edu]
Sent: Thursday, September 30, 2010 3:09 PM
To: Ravi Varadhan
Cc: arindam fadikar; r-help@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: rvarad...@jhmi.edu
----- Original Message -----
From: Ravi Varadhan <rvarad...@jhmi.edu>
Date: Thursday, September 30, 2010 3:04 pm
Subject: Re: [R] how to avoid NaN in optim()
To: Ravi Varadhan <rvarad...@jhmi.edu>
Cc: arindam fadikar <arindam.fadi...@gmail.com>, r-help@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: rvarad...@jhmi.edu
>
>
> ----- Original Message -----
> From: Ravi Varadhan <rvarad...@jhmi.edu>
> Date: Thursday, September 30, 2010 2:54 pm
> Subject: Re: [R] how to avoid NaN in optim()
> To: arindam fadikar <arindam.fadi...@gmail.com>
> Cc: r-help@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: rvarad...@jhmi.edu
> >
> >
> > ----- Original Message -----
> > From: arindam fadikar <arindam.fadi...@gmail.com>
> > Date: Thursday, September 30, 2010 2:17 pm
> > Subject: [R] how to avoid NaN in optim()
> > To: r-help@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@r-project.org mailing list
> > >
> > > PLEASE do read the posting guide
> > > and provide commented, minimal, self-contained, reproducible code.
> >
> > ______________________________________________
> > R-help@r-project.org mailing list
> >
> > PLEASE do read the posting guide
> > and provide commented, minimal, self-contained, reproducible code.
>
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
