Guest User

Untitled

a guest
Jan 19th, 2018
88
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 0.92 KB | None | 0 0
  1. library(rootSolve)
  2.  
  3. d <- read.table(text="indx rate n d
  4. 1 0.12 158 14
  5. 2 0.095 135 9
  6. 3 0.057 123 4
  7. 4 0.033 115 5
  8. 5 0.019 90 4", header=T)
  9.  
  10. d$real <- with(d, d/n)
  11.  
  12. opt <- d[ ,c("rate","real", "n")]
  13.  
  14. # this is close to the correct solution!
  15.  
  16. scaler <- apply(opt, 1, function(z) uniroot.all(
  17. function(alpha) z[2] - (1 / (1 + alpha * ( (1 - z[1]) / z[1] )) ), interval = c(0,10)))
  18.  
  19. # check for solution (not fully correct!)
  20.  
  21. round(crossprod(scaler * opt$real, opt$n)/sum(opt$n), 3) == round(crossprod(round(opt$rate, 3), opt$n)/sum(opt$n), 3)
  22.  
  23.  
  24. # using optim() - completely wrong results
  25.  
  26. infun <- function(data, alpha){ l <- with(data,
  27. ( rate - (1 / ( 1 + alpha[1] * ( (1 - real)/real ))) )); return( -sum( l ) ) }
  28.  
  29. opt_out <- optim(c(0,0), infun, data=opt, method = "BFGS", hessian = TRUE)
  30.  
  31. with(opt, (1 / ( 1 + opt_out$par[1] * ( (1 - real)/real ))))
Add Comment
Please, Sign In to add comment