CODEWITHSUNDEEP
roptimization
Maximilian
Maximilian

Reputation: 4229

How to find solution with optim using r

I have quite simple question yet I cannot seem to solve it.

I'm trying to find a parameter with optim() r function.

Here is the case:

library(rootSolve)

d <- read.table(text="indx rate n d 
            1   0.12    158 14
            2   0.095   135 9
            3  0.057    123 4
            4  0.033    115 5
            5    0.019  90  4", header=T)

d$real <- with(d, d/n)

opt <- d[ ,c("rate","real", "n")]

# this is close to the correct solution!

scaler <- apply(opt, 1, function(z) uniroot.all( 
function(alpha) z[2] - (1 / (1 + alpha * ( (1 - z[1]) / z[1] )) ), interval = c(0,10)))

# check for solution (not fully correct!)

round(crossprod(scaler * opt$real, opt$n)/sum(opt$n), 3) == round(crossprod(round(opt$rate, 3), opt$n)/sum(opt$n), 3)


# using optim() - completely wrong results

infun <- function(data, alpha){ l <- with(data, 
( rate - (1 / ( 1 + alpha[1] * ( (1 - real)/real ))) )); return( -sum( l ) ) }

opt_out <- optim(c(0,0), infun, data=opt, method = "BFGS", hessian = TRUE)

with(opt, (1 / ( 1 + opt_out$par[1] * ( (1 - real)/real ))))

Upvotes: 0

Views: 467

Answers (1)

MrSmithGoesToWashington
MrSmithGoesToWashington

Reputation: 1076

You are trying, with your code, to get an unique alpha for all, but you want to have five values .. So, you are leaded to make a sum .. but if you have negative and positive values, your sum could go near zero even with individual terms far from 0 ..

Moreover, your infun function is not in accordance with your previous function ..

What you can do is something like that : 

infun <- function(alpha){ l <- with(cbind(d, alpha), ( real - (1 / ( 1 + 
alpha * ( (1 - rate)/rate ))) )); return( sum(abs(l)) ) }
param <-  c(5,5,5,5,5)
opt_out <- optim(par = scaler, infun,   method = "BFGS", hessian = TRUE)

And in order to check the result you should have written : 

    with( cbind(opt,opt_out$par), real -1 / ( 1 + opt_out$par * ( (1 - rate)/rate )))

To get the true solution, you can do (after a very litle mathematics on a paper) : 

sol <- -((opt[,2]-1)/(opt[,2]))*(opt[,1]/(1-opt[,1]))

and test it :

with( cbind(opt,sol), real -1  / ( 1 + opt_out$par * ( (1 - rate)/rate )))

Upvotes: 1

Related Questions