Kate,
This is a difficult problem, mainly because your function is not
deterministic. Run the function a couple of times with the same parameter
values and you will get a different value each time. Obviously, you cannot
optimize such functions in the ususal sense.
However, you can rewrite the function as follows:
mu2 <- 0.4
sigma2 <- 4
tau <- 0.75
mean.diff <- 2
var.ratio <- 3
set.seed(128)
u <- runif(1e06)
Y <- qnorm(u,mean=mu2,sd=sigma2)
mY <- mean(Y)
sY <- sd(Y)
u <- runif(1e06)
X0 <- qnorm(u,mean=mu2,sd=sigma2)
parameter2 <- function(cons) {
a<-cons[1]
b<-cons[2]
f <- rep(NA, 2)
X <- X0 + a*abs(u-tau)^b*(u>tau)
f[1] <- abs(mean(X) - mY) - mean.diff
f[2] <- sd(X) - sqrt(var.ratio) * sY
f
}
Now, parameter2() is a deterministic function because it is conditioned upon
X and Y, i.e. X and Y are fixed now. Still this is not a trivial problem.
The equations are non-smooth, so most standard optimization techniques will
struggle with it. However, I have had success with using the dfsane()
function in "BB" package to solve non-smooth systems arising in rank-based
regression methods in survival analysis. It also seems to work for your
example:
require(BB) # you have to install BB package
c0 <- c(3,-1)
ans <- dfsane(par=c0, fn=parameter2, method=3, control=list(NM=TRUE)) #
this converges to a solution
> ans
$par
[1] 4.9158558 -0.1941085
$residual
[1] 8.17635e-10
$fn.reduction
[1] 0.3131852
$feval
[1] 301
$iter
[1] 158
$convergence
[1] 0
$message
[1] "Successful convergence"
The solution, of course, depends upon X and Y. So, if you generate another
set of X and Y, you will get a different answer.
Hope this helps,
Ravi.
----------------------------------------------------------------------------
-------
Ravi Varadhan, Ph.D.
Assistant Professor, The Center on Aging and Health
Division of Geriatric Medicine and Gerontology
Johns Hopkins University
Ph: (410) 502-2619
Fax: (410) 614-9625
Email: rvarad...@jhmi.edu
Webpage:
http://www.jhsph.edu/agingandhealth/People/Faculty_personal_pages/Varadhan.h
tml
----------------------------------------------------------------------------
--------
-----Original Message-----
From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On
Behalf Of yh...@illinois.edu
Sent: Friday, July 17, 2009 10:25 AM
To: r-help@r-project.org
Subject: Re: [R] Solving two nonlinear equations with two knowns
My question is as follows
Y~N(mu2,sigma2^2), i.e. Y has cdf F2
X generates from:
X=F2^{-1}(u)+a*|u-tau|^b*I(u>0.75), where u~U(0,1)
Given tau=0.75, I want to find a and b such that
E(X)-E(Y)=2 and Var(X)/Var(Y)=3
I try optim and grid seaching and get different results, any solution?
Thanks,
Kate
##########
#R code:
mu2=0.4
sigma2=4
tau=0.75
mean.diff=2
var.ratio=3
#Use optim:
parameter<-function(c)
{
a<-c[1]
b<-c[2]
u<-runif(10000)
Y<-qnorm(u,mean=mu2,sd=sigma2)
u<-runif(10000)
X<-qnorm(u,mean=mu2,sd=sigma2)+a*abs(u-tau)^b*(u>tau)
return((abs(mean(X)-mean(Y))-mean.diff)^2+(var(X)/var(Y)-var.ratio)^2)
}
c0<-c(3,1)
cstar<-optim(c0,parameter)$par
astar<-cstar[1] #4.1709
bstar<-cstar[2] #-0.2578
#Use grid seaching (I randomly assign a rage (-10,100)):
parameter.X<-function(a,b)
{
TSE<-matrix(0, length(a),length(b))
u<-runif(10000)
Y<-qnorm(u,mean=mu2,sd=sigma2)
u<-runif(10000)
for(i in 1: length(a))
{
for(j in 1: length(b))
{
X<-qnorm(u,mean=mu2,sd=sigma2)+a[i]*abs(u-tau)^b[j]*(u>tau)
TSE[i,j]<-(abs(mean(X)-mean(Y))-mean.diff)^2+(var(X)/var(Y)-var.ratio)^2
}
}
minTSE<-min(TSE)
a.optimal<-a[which(TSE==min(TSE),arr.ind = TRUE)[1]]
b.optimal<-b[which(TSE==min(TSE),arr.ind = TRUE)[2]]
return(list(a.optimal=a.optimal,b.optimal=b.optimal,minTSE=minTSE))
}
a0<-seq(-10,100,,50)
b0<-seq(-10,100,,50)
tse1<-parameter.X(a0,b0)
astar<-tse1$a.optimal # 84.28571
bstar<-tse1$b.optimal #1.224490
______________________________________________
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.
______________________________________________
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.
