[R] Solving two nonlinear equations with two knowns

[Ravi Varadhan]

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: rvaradhan at jhmi.edu

Webpage:
http://www.jhsph.edu/agingandhealth/People/Faculty_personal_pages/Varadhan.h
tml

 

----------------------------------------------------------------------------
--------


-----Original Message-----
From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org] On
Behalf Of yhsu6 at illinois.edu
Sent: Friday, July 17, 2009 10:25 AM
To: r-help at 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 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.