[R] using optimize with two unknowns, e.g. to parameterize a distribution with given confidence interval

[Ravi Varadhan]

You cannot use `optimize' when there are two or more parameters to be
optimized.  I don’t know if other have suggested any solution to this, but
here are 2 approaches:

# Estimating LCL and UCL separately using `optimize'.

prior.lcl <- function(x, alpha, mean, var) {
  a <- abs(plnorm(x, mean, var) - (alpha/2))
  return(a)
  }

prior.ucl <- function(x, alpha, mean, var) {
  b <- abs(plnorm(x, mean, var) - (1-alpha/2))
  return(b)
  }


optimize(f=prior.lcl, mean=0, var=1, alpha=0.05, interval=c(0,20))

optimize(f=prior.ucl, mean=0, var=1, alpha=0.05, interval=c(0,20))


# Combining LCL and UCL estimation
# Using `optimx' package to illustrate how to solve this using different
optimizers
# This also shows that this problem is not so easy to solve
#
prior <- function(x, alpha, mean, var) {
  a <- abs(plnorm(x[1], mean, var) - (alpha/2))
  b <- abs(plnorm(x[2], mean, var) - (1-alpha/2))
  return(a+b)
  }


require(optimx)

optimx(par=c(0,10), fn=prior, mean=0, var=1, alpha=0.05, method=c("Nelder",
"BFGS", "CG", "spg", "bobyqa", "Rvmmin", "nlminb", "Rcgmin", "ucminf"))

optimx(par=c(1,10), fn=prior, mean=0, var=1, alpha=0.05, method=c("Nelder",
"BFGS", "CG", "spg", "bobyqa", "Rvmmin", "nlminb", "Rcgmin", "ucminf"))


Hope this helps,
Ravi.


-----Original Message-----
From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org] On
Behalf Of David LeBauer
Sent: Friday, October 15, 2010 5:31 PM
To: r-help
Subject: [R] using optimize with two unknowns, e.g. to parameterize a
distribution with given confidence interval

Hi,

I would like to write a function that finds parameters of a log-normal
distribution with a 1-alpha CI of (x_lcl, x_ucl):

However, I don't know how to optimize for the two unknown parameters.

Here is my unsuccessful attempt to find a lognormal distribution with
a 90%CI of 1,20:

prior <- function(x_lcl, x_ucl, alpha, mean, var) {
  a <- (plnorm(x_lcl, mean, var) - (alpha/2))^2
  b <- (plnorm(x_ucl, mean, var) - (1-alpha/2))^2
  return(a+b)
  }

optimize(fn=prior, interval = c(-5, 100), 1, 20)

I understand that this problem has a closed form solution, but I would
like to make this a general function.

Thanks,

David

-- 
David LeBauer, PhD
Energy Biosciences Institute
University of Illinois Urbana-Champaign
1206 W. Gregory Drive
Urbana, IL  61801, U.S.A.

______________________________________________
R-help at 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.