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.
______________________________________________
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.
