# MLE for parametric binomial model

I have a model in which $$p_i=f(\theta,Z_i)$$, where $$Z_i$$ are iid latent variables distributed with CDF $$F_\theta$$, and $$d_i\sim B(n_i,p_i)$$, where $$B$$ is the binomial distribution. The likelihood function is thus (up to a constant) $$L(\theta\mid d_1,\dots,d_k)=\int_{\mathbb{R}}\dots \int_{\mathbb{R}}\prod_{i=1}^k f(\theta,z_i)^{d_i}(1-f(\theta,z_i))^{n_i-d_i} dF_\theta(z_1)\dots dF_\theta(z_k).$$ The problem is that when $$p_i$$ fluctuates around $$0.05$$, and the $$n_i$$ are around $$1000$$ each, and $$k$$ is around $$20$$, then the magnitude of the likelihood is so small that numerically evaluating these integral becomes impossible (it seems). Is there some way around this?

• Depending how big $k$ you are likely suffering the curse of dimensionality. The issue then has little to do with the "fluctuating" value of $p_i$ but rather that you're integrating over a lot of stuff. It looks a good application for INLA if you're willing to go Bayesian (which is an ideal framework for dealing with latent variables). Commented Jul 20, 2023 at 20:53
• @AdamO Thanks, not familiar with INLA but it looks really interesting. I'll see if it can help me further. Commented Jul 21, 2023 at 14:02
• Are you using software where you can access and control multiple-precision arithmetic and integration routines?
– JimB
Commented Aug 3, 2023 at 4:27
• Can you not switch the order of the product and integration? Then you can look at the log of the likelihood?
– JimB
Commented Aug 3, 2023 at 4:52

Is the following a specific example of what you're trying to do? If so, dealing with the log of the likelihood is the way to avoid over and/or underflow issues.

Suppose $$Z_i$$ has a beta distribution ($$Z_i \sim B(a,b)$$) and $$f(\theta,Z_i)=Z_i$$.

$$L(a,b|d_i)=\int_0^1 z^{d_i}(1-z)^{n_i-d_i} \frac{z^{a-1} (1-z)^{b-1} \Gamma (a+b)}{ \Gamma (a) \Gamma (b)} dz=\frac{\Gamma (a+b) \Gamma (a+d_i) \Gamma (n_i-d_i+b)}{\Gamma (a) \Gamma (b) \Gamma (n_i+a+b)}$$

$$\log{L(a,b|d_1,\ldots,d_k)}=\sum_{i=1}^k \log\left(\frac{\Gamma (a+b) \Gamma (a+d_i) \Gamma (n_i-d_i+b)}{\Gamma (a) \Gamma (b) \Gamma (n_i+a+b)}\right)$$

Using Mathematica the maximum likelihood estimates can be obtained with the following:

n = 1000;
k = 20;
SeedRandom[12345];
d = RandomVariate[BinomialDistribution[n, #], 1][[1]] & /@ p;
logL = Sum[LogGamma[a + b] + LogGamma[a + d[[i]]] + LogGamma[n + b - d[[i]]] -
LogGamma[a] - LogGamma[b] - LogGamma[n + a + b], {i, 1, k}];
FindMaximum[{logL, a > 0 && b > 0}, {{a, 3}, {b, 60}}]
(* {-4116.13, {a -> 3.20512, b -> 55.2461}} *)


Using R:

# Set sample sizes and values for Z and p and
# get random samples
set.seed(12345)
k <- 20
n <- rpois(k, 1000)
z <- rbeta(k, 3, 57)
d <- NULL
for (i in 1:k) {
d[i] <- rbinom(1, n[i], z[i])
}

# Log of the likelihood function
logL <- function(parms, k=k, n=n, d=d) {
a <- parms[1]
b <- parms[2]
sum(lgamma(a + b) - lgamma(a) - lgamma(b) +
lgamma(d + a) + lgamma(n - d + b) - lgamma(n + a + b))
}

# Maximum likelihood estimates
mle <- optim(c(3,60), logL, k=k, n=n, d=d, control=list(fnscale=-1))
mle
$par [1] 2.704537 47.144652$value
[1] -4118.709

$counts function gradient 55 NA$convergence
[1] 0

\$message
NULL