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]]
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.
I forgot to mention: You actually got correct results with using optim and `CG'. The warning messages were just telling you that there were some log(negative number) operations during the iterative process. 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.
arindam fadikar wrote:> > > 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 > } > ..... >Extending the tests in the if in loglik to if (c(p,q,r) > rep(0,3) && c(p,q,r) < rep(1,3) && (p^2 + 2*p*r)>0 && (q^2 + 2*q*r)>0) would also help. /Berend -- View this message in context: http://r.789695.n4.nabble.com/how-to-avoid-NaN-in-optim-tp2738093p2746635.html Sent from the R help mailing list archive at Nabble.com.