CODEWITHSUNDEEP
roptimization
MLEN
MLEN

Reputation: 2561

Using optimisation function (optimise) together with dbinom in R (optimisation issue)

When p = 0.5, n = 5 and x = 3

dbinom(3,5,0.5) = 0.3125

Lets say I dont know p (n and x is known) and want to find it. 

binp <- function(bp) dbinom(3,5,bp) - 0.3125
optimise(binp, c(0,1))

It does not return 0.5. Also, why is

dbinom(3,5,0.5) == 0.3125 #FALSE

But,

x <- dbinom(3,5,0.5)
x == dbinom(3,5,0.5) #TRUE

Upvotes: 0

Views: 128

Answers (1)

cuttlefish44
cuttlefish44

Reputation: 6796

optimize() searches the param that minimizes the output of function. Your function can return a negative value (e.g., binp(0.1) is -0.3044). If you search the param that minimizes difference from zero, it would be good idea to use sqrt((...)^2). If you want the param that makes output zero, uniroot would help you. And the param what you want isn't uniquely decided. (note; x <- dbinom(3, 5, 0.5); x == dbinom(3, 5, 0.5) is equibalent to dbinom(3, 5, 0.5) == dbinom(3, 5, 0.5))

 ## check output of dbinom(3, 5, prob)
input <- seq(0, 1, 0.001)
output <- Vectorize(dbinom, "prob")(3, 5, input)
plot(input, output, type="l")
abline(h = dbinom(3, 5, 0.5), col = 2)     # there are two answers

enter image description here

max <- optimize(function(x) dbinom(3, 5, x), c(0, 1), maximum = T)$maximum  # [1] 0.6000006

binp <- function(bp) dbinom(3,5,bp) - 0.3125   # your function

uniroot(binp, c(0, max))$root # [1] 0.5000036
uniroot(binp, c(max, 1))$root # [1] 0.6946854


binp2 <- function(bp) sqrt((dbinom(3,5,bp) - 0.3125)^2)

optimize(binp2, c(0, max))$minimum  # [1] 0.499986
optimize(binp2, c(max, 1))$minimum  # [1] 0.6947186


dbinom(3, 5, 0.5) == 0.3125            # [1] FALSE
round(dbinom(3, 5, 0.5), 4) == 0.3125  # [1] TRUE
format(dbinom(3, 5, 0.5), digits = 16) # [1] "0.3124999999999999"

Upvotes: 2

Related Questions