Calculating an integral of a Bernoulli likelihood function in R

I am trying to integrate a very simple likelihood function (three Bernoulli trials) in R, but it seems that I have not understood how integrate function works. My code is as follows:

integrate(
    function(theta) prod(
        dbinom(
            x = c(1,0,1), #x is the data: "success, fail, success"
            size = 1, #Just one Bernoulli trial times the length of x.
            prob = theta, #Since, this is a likelihood function, the parameter value is the variable
            log = F #We use a likelihood instead of log-likelihood
        )
    ),
    lower = 0, #The parameter theta cannot be smaller than zero...
    upper = 1 #...nor greater than one.
)

However, I get this error message:

Error in integrate(function(theta) prod(dbinom(x = c(1, 0, 1), size = 1,  : 
  evaluation of function gave a result of wrong length

Now, why is the result of wrong length? The result is always of length 1, since the anonymous function uses prod function, which in turn creates a product of all the vector elements the function dbinom returns (this vector has the length of three, since its first argument x has length of three).

What the result should be then to be of right length?

Solution

The issue here is not with dbinom but with prod. dbinom is a vectorized function but prod, according to its definition, returns a "a numeric (of type "double") or complex vector of length one."

As noted in the comments, the simplest approach is to simply wrap your function in Vectorize inside the integrate call. For example:

fn <- function(theta) prod(dbinom(c(1, 0, 1), size = 1, prob = theta, log = FALSE))
integrate(Vectorize(fn), 0, 1)
0.08333333 with absolute error < 9.3e-16