2
$\begingroup$

Suppose there is a hidden gaussian with mean $\mu$ and variance $\sigma^2$, and that $X_i \sim \mathcal{N}(\mu,\sigma^2)$ where the $X_i$ are i.i.d. If I can only oberve the rounded value of $X_i$, i.e. $Y_i = \lfloor X_i + 1/2 \rfloor$, is there an effective means of computing the maximum likelihood estimates of $\mu$ and $\sigma$?

$\endgroup$
  • 1
    $\begingroup$ Maximize the log of the likelihood: $\log L=\sum _{i=1}^n \log \left(\Phi \left(\frac{\left(y_i-\mu +\frac{1}{2}\right){}^2}{2 \sigma ^2}\right)-\Phi \left(\frac{\left(y_i-\mu -\frac{1}{2}\right){}^2}{2 \sigma ^2}\right)\right)$. $\endgroup$ – JimB Feb 5 at 5:18
  • $\begingroup$ Any statistical software they let's you specify that data are censored to lie in intervals or that lets you specify your own likelihood should do the job. $\endgroup$ – Björn Feb 5 at 5:34
  • $\begingroup$ That makes sense. I have accepted the answer from @JimB. I am in a situation where I want to update mu and sigma for 50000 different gaussians 15 times a second on commodity hardware, so I doubt that a numerical solver will be fast enough. $\endgroup$ – Josh Brown Kramer Feb 7 at 16:36
  • $\begingroup$ Not sure what you mean by "I doubt that a numerical solver will be fast enough." There is no simple closed-form solution for such censored data. All approaches will be iterative. But if there's a wide range of integers, then just calculating the usual sample mean and standard deviation might be good enough. (And "good enough" is a subject matter issue and not a statistical issue.) $\endgroup$ – JimB Feb 7 at 17:19
1
$\begingroup$

As @Björn suggests, using package software is likely to work best in that the code is (again, more likely to be) debugged, have better documentation, be repeatable, less prone to numerical precision issues, and gives you someone else to blame.

But if you need to write your here is a bare-bones approach using R:

# Generate data
  set.seed(12345)
  x = rnorm(100, mean=20, sd=2)

# Round data
  y = round(x)

# Define log likelihood function
  logL = function(p, z) {
    mu = p[1]
    sigma = p[2]
    sum(log(pnorm(z+0.5, mu, sigma)-pnorm(z-0.5, mu, sigma)))
  }

# Initial values
  mu0 = mean(y)
  sigma0 = sd(y)

# Find maximum likelihood estimates
  sol = optim(c(mu0, sigma0), logL, z=y, control=list(fnscale=-1), hessian=TRUE)

# Show results
  cat("MLE of mu: ", sol$par[1], "\n")
# MLE of mu:  20.46009

  cat("MLE of sigma:", sol$par[2], "\n")
# MLE of sigma: 2.277309 

# Estimates of standard errors and covariance
  cov = solve(-sol$hessian)
  cat("se of MLE of mu:", cov[1,1]^0.5, "\n")
# se of MLE of mu: 0.2295532 
  cat("se of MLE of sigma:", cov[2,2]^0.5, "\n")
# se of MLE of sigma: 0.1636504 
$\endgroup$
  • $\begingroup$ Note that this can be made more efficient if the frequencies of each integer are tabulated. But unless you're also doing a bootstrap or have a huge sample size, there might not be much improvement in speed. $\endgroup$ – JimB Feb 5 at 16:23
  • $\begingroup$ Nice touch with the Hessian. I think that is in the spirit of what the OP was asking. :) $\endgroup$ – Peter Leopold Feb 5 at 17:20
  • $\begingroup$ @PeterLeopold. Gotta have measures of precision with every estimate. One of my standard soapbox speeches/sermons. $\endgroup$ – JimB Feb 5 at 18:03
0
$\begingroup$

(My solution is basically the same as what @JimB suggested in his comment.)

First, some theory:

We can think of $[x]=x+\epsilon$ where $f(\epsilon) \sim (1-\delta)+2\delta\epsilon$ for $1/2\le \epsilon \le 1/2$. This is a sloping distribution that approximates $N(x|\mu,\sigma^2)$) by a line with slope $2\delta=N'([x]|\mu,\sigma^2)=-\frac{[x]-\mu}{\sigma^2} N([x]|\mu,\sigma^2)$. We see that $\delta$ is nearly zero when $x \sim \mu$ or when $x\gg\mu$ or $x\ll\mu$. $\delta$ is a maximum value around the inflection points $[x] \sim \pm \sigma$, which is where a lot of the action happens.

The mean of $y$ is found to be an unbiased estimator of $\mu$. $$E(y)=E([x])=E(x+\epsilon)=E(x)+E(\epsilon)=E(x)=\mu$$

where I have used $E(\epsilon)=0$ since $E(\delta)=0$ by symmetry.

For variance,

$$var(y)=var([x])=var(x+\epsilon)=var(x)+var(\epsilon)$$

where $var(\epsilon)$ is estimated as a function of $\delta$ to be

$$var(\epsilon) \sim 1/12 - \delta^2/36=\frac{1}{12} (1-\delta^2/3).$$

This would now have to be averaged over the distribution of values of $\delta,$ which is doable but tedious.

On 10,000 simulations of 1000 random normals with means around 50 and variances around 100, I find

enter image description here

which contains at least one big surprise. The average value is $1/12$ as perhaps expected, but this figure shows that for about half of the simulations, the effect of rounding the data to integer is to underestimate $\sigma^2$ rather than to overestimate it, as would be expected from basic probability theory.

Perhaps this is because there is assumption that $x$ and $\epsilon$ are independent is false!

Regardless, now to the optimization problem in $\mu$ and $\sigma$

Here is an example, in R.

rm(list=ls())

mu=50
sigma=10
N=1000
data=rnorm(N,mu,sigma)
idata=round(data,0)

indices<-seq(round(mean(idata)-5*sd(idata)),round(mean(idata)+5*sd(idata)))

logLike<-function(pars,idata,indices){
  mu=pars[1]
  sigma=pars[2]
  imin=min(indices)
  ps<-pnorm(indices+.5,mu,sigma)-pnorm(indices-.5,mu,sigma)
  logps<-log(ps)
  -sum(logps[idata-imin+1])
}
startingPars=c(mean(idata),sd(idata))
oo<-optim(startingPars,logLike,idata=idata,indices=indices,control=list(reltol=1e-12))
print(paste("target Sigma",sigma))
print(paste("sample sigma",sd(data)))
print(paste("sample integerized sigma",sd(idata)))
print(paste("optimized sigma",oo$par[2]))

One example of this gives

[1] "target variance 127.743486150357"
[1] "sample variance 124.050532385344"
[1] "sample rounded variance 124.173364364364"
[1] "optimized variance 123.965644616766"

In this one case, we see that rounding increases the variance and that optimizing reduces the estimate to a value closer to the original sample variance, which is the best linear unbiased estimator of the target variance.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.