# [R] (no subject)

Berend Hasselman bhh at xs4all.nl
Mon Sep 12 08:54:00 CEST 2011

```li li-13 wrote:
>
> Dear all,
>     Can anyone take a look at my program below?
> There are two functions: f1 (lambda,z,p1) and f2(p1,cl, cu).
> I fixed p1=0.15 for both functions. For any fixed value of lambda (between
> 0.01 and 0.99),
> I solve f1(p1=0.15, lambda=lambda, z)=0 for the corresponding cl and cu
> values.
> Then I plug the calculated cl and cu back into the function f2.
> Eventually, I want to find the lambda value and the corresponding cl and
> cu
> values that would
> make f2=0.1.
>    The result of this program does not seem to match the answer I have.
> Can
> some one give me
> some hint? Thank you very much.
>             Hannah
>
>
>
> u1 <- -3
>
> u2 <- 4
>
>
> f1 <- function(lambda,z,p1){
>
> lambda*(p1*exp(u1*z-u1^2/2)+(0.2-p1)*exp(u2*z-u2^2/2))-(1-lambda)*0.8}
>
>
> f2 <- function(p1,cl, cu){
>
>
> 0.8*(pnorm(cl)+(1-pnorm(cu)))/(0.8*(pnorm(cl)+(1-pnorm(cu)))+p1*(pnorm(cl-
> u1)+(1-pnorm(cu-u1)))+(0.2-p1)*(pnorm(cl-u2)+(1-pnorm(cu-u2))))}
>
>
> p1 <- 0.15
>
>
> lam <- seq(0.01,0.99, by=0.001)
>
> x1 <- numeric(length(lam))
>
>
> for (i in 1:length(lam)){
>
>
>
> cl <- uniroot(f1, lower =-10, upper = 0,
>
>            tol = 1e-10,p1=p1,lambda=lam[i])\$root
>
>
> cu <- uniroot(f1, lower =0, upper = 10,
>
>            tol = 1e-10,p1=p1,lambda=lam[i])\$root
>
>
>
> x1[i]<- f2(p1=p1, cl=cl, cu=cu) }
>
>
>
> k <- 1
>
> while(k<length(lam) && x1[k]<=0.1){
>
>     k=k+1
>
>   }
>
>   k<-k-1;k
>
> lower <- uniroot(f1, lower =-10, upper = 0,
>
>            tol = 1e-10,p1=p1,lambda=lam[k])\$root
>
> upper <- uniroot(f1, lower =0, upper = 10,
>
>            tol = 1e-10,p1=p1,lambda=lam[k])\$root
>
> res <- c(lower, upper)
>

properly formatted code and a clear description of what you expect and what
you are getting.

Since I can't resist the temptation.

Method 1:
-----------

u1 <- -3
u2 <- 4

f1 <- function(lambda,z,p1){
lambda*(p1*exp(u1*z-u1^2/2)+(0.2-p1)*exp(u2*z-u2^2/2))-(1-lambda)*0.8 }

f2 <- function(p1,cl, cu){

0.8*(pnorm(cl)+(1-pnorm(cu)))/(0.8*(pnorm(cl)+(1-pnorm(cu)))+p1*(pnorm(cl-
u1)+
(1-pnorm(cu-u1)))+(0.2-p1)*(pnorm(cl-u2)+(1-pnorm(cu-u2))))
}

p1 <- 0.15
lam <- seq(0.01,0.99, by=0.001)

df.c <- data.frame(cl=numeric(length(lam)),cu=numeric(length(lam)))

for (i in 1:length(lam)){
cl <- uniroot(f1, lower =-10, upper = 0, tol =
1e-10,p1=p1,lambda=lam[i])\$root
cu <- uniroot(f1, lower =0, upper = 10,  tol =
1e-10,p1=p1,lambda=lam[i])\$root
df.c[i,"cl"] <- cl
df.c[i,"cu"] <- cu
}

x1 <- f2(p1=p1,cl=df.c[,"cl"], cu=df.c[,"cu"])
df.c[,"x1"] <- x1

df.c

# find the index for which the deviation from the target value 0.1 is
smallest
x.minidx <- which.min(abs(x1-.1))
x.minidx
x1[x.minidx]

Method 2:
------------
But why so complicated? You want a lambda that makes f2 as close as possible
to .1.

Define a target function  with one argument lambda

target <- function(lambda) {
cl <- uniroot(f1, lower =-10, upper = 0, tol =
1e-10,p1=p1,lambda=lambda)\$root
cu <- uniroot(f1, lower =0, upper = 10,  tol =
1e-10,p1=p1,lambda=lambda)\$root
f2(p1=p1, cl=cl, cu=cu) - 0.1
}

and solve for lambda as follows

res <- uniroot(target, lower=0.01, upper=0.99)
res
lambda <- res\$root

# Show some results

cl <- uniroot(f1, lower =-10, upper = 0, tol =
1e-10,p1=p1,lambda=lambda)\$root
cu <- uniroot(f1, lower =0, upper = 10,  tol =
1e-10,p1=p1,lambda=lambda)\$root
f2(p1=p1, cl=cl, cu=cu)

This is more accurate than Method 1.

/Berend

--
View this message in context: http://r.789695.n4.nabble.com/no-subject-tp3805801p3806668.html
Sent from the R help mailing list archive at Nabble.com.

```