[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)
> 

You should follow the advice of the previous repliers: meaningful subject,
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.



More information about the R-help mailing list