# Estimating the mean of a truncated gaussian curve

Say I have a black box generating data samples, and I want to estimate the parameters of the black box from the samples.

The black box works like this: it has a parameter m (a real number), and to generate a value v, it first generates v0 according to a normal distribution (with mean m and variance 1), and if v0 is positive it returns v0, if not it returns 0.

So my data samples will be a bunch of zeroes and positive real numbers.

So my question is, from a sample, how do I estimate m?

And what kind of mathematical tools do I use to reason about this?

To me, this looks like a straightforward case of bayesian probability, where I would use p(samples|m) to get p(m|samples) and have some prior on the distribution of m.

So uncle Bayes would say: $p(m|series) = p(series|m) * p(m) / p(series)$

Since the samples are independant, $p(samples|m) = \prod p(sample|m)$

...but some of those $p(sample|m)$ are "discrete probabilities" (when the value is 0), and some are continuous probabilities! Can I muliply them just like that?

(Same goes for calculating $p(samples)$)

Can someone help me clear the confusion?

-
I assume m is the variance of the normal distribution? – Vectornaut Apr 14 '10 at 16:01
No, the mean actually, sorry :P (edited now) – Emile Apr 14 '10 at 16:51

OK, let me fully address the question since there is no easy way out. The normal approach is to maximize the "likelihood" of the data under the parameter. The key question here is how to define likelihood for a mixed distribution. Let's use the standard approach as our guide.

Parameter estimation is usually based on the idea that we want to choose parameters that make our data "the most likely." For a discrete probability distribution, we interpret this to mean that our data is the most probable. But this breaks down in the case of continuous probability distributions, where, no matter our choice of parameters, our data has probability zero.

Statisticians thus replace the probability with the probability density for continuous distributions. Here is the justification for this. Instead of actually having a set of numbers drawn from the probability distribution, you have a highly accurate measurement---say, your sequence $\{x_i\}$ for $i = 1,\dots,n$ tells you that the true value of the (still unknown) sequence $\{g_i\}$ satisfies $|x_i - g_i| < \varepsilon$ for all $i$. When $\varepsilon$ is sufficiently small, the replacement $$\mathbb{P}(|x_i - g_i|) < \varepsilon )\approx \varepsilon p_{g}(x_i)$$ is very accurate, where $p_g$ is the pdf of $g_i$. Assuming that your sequence is iid, we are led to the approximation $$\mathbb{P}(|x_i - g_i| < \varepsilon \text{ for all } i) \approx \varepsilon^n \prod_{i=1}^n p_g(x_i).$$ We thus choose the pdf from our family which maximizes the right hand side of the above equation, reproducing the standard maximum likelihood method.

Now the question is, what do we do with mixed distributions? When there is a mass at a point $x_i$, that is $\mathbb{P}(x_i=g_i) > 0$, our first approximation is incorrect; for very small $\varepsilon$, we have the approximation $$\mathbb{P}(|x_i - g_i| < \varepsilon) \approx \mathbb{P}(x_i = g_i)$$ If we let $\mathcal{N}$ be the index set of the "massless" samples, we can approximate the probability of our data as $$\mathbb{P}(|x_i - g_i| < \varepsilon) \approx \varepsilon^n \prod_{i \in \mathcal{N}} p_g(x_i) \prod_{i \notin \mathcal{N}} \mathbb{P}(x_i = g_i).$$ where $n$ is the number of elements in $\mathcal{N}$. That is, we can reasonably define our maximum likelihood estimate for a parameter $m$ as the value of the parameter that maximizes $$\prod_{i \in \mathcal{N}} p_g(x_i) \prod_{i \notin \mathcal{N}} \mathbb{P}(x_i = g_i).$$

In your case, it is fairly simple to write down the value of the likelihood function above. First, note that $$\mathbb{P}(x=0) = \frac{1}{\sqrt{2\pi}} \int_{-\infty}^{-m} e^{-x^2/2}dx.$$ For $x>0$, you have the standard Gaussian pdf $p_g(x) = \tfrac{1}{\sqrt{2\pi}} e^{-(x-m)^2/2}$.

I won't do any more here; suffice it to say that the standard approach to maximizing the likelihood involves taking the logarithm of the likelihood function and setting its derivative to zero. You will probably get a transcendental equation that you will need to solve numerically.

-

I posted a a follow up question here and http://stats.stackexchange.com/questions/2037/estimating-mean-and-st-dev-of-a-truncated-gaussian-curve-without-spike; I got this answer in return and thought it may be found interesting: Cohen and Hald solved this problem using ML (with a Newton-Raphson root finder) around 1950 - http://projecteuclid.org/euclid.aoms/1177729751 (open access). Another paper is Max Halperin's "Estimation in the Truncated Normal Distribution" available on JSTOR http://www.jstor.org/stable/2281315 (for those with access). Googling "truncated gaussian estimation" produces lots of useful-looking hits.

-

Edit: This is an answer to a different question (see Emile's comment).

I don't think the question you asked gets at the issue you want to talk about, because, as Mike McCoy said, the zero samples are irrelevant to the parameter estimation. You can just ignore the zeroes and estimate $m$ from the nonzero samples.

Here's an alternative question, which might work better. Say you have a white box that works like this. It has a parameter $x^\ast$, and to generate a value, it first flips a fair coin. If the coin comes up heads, the box returns $x^\ast$. If the coin comes up tails, the box returns a sample from the standard normal distribution. How do we estimate $x^\ast$?

This problem has a strange feature. If we have the ability to compare two real numbers exactly, we can find $x^\ast$ exactly from just four samples, on average! The probability of getting two equal samples from a normal distribution is zero; so if we get two equal samples, they're equal to $x^\ast$ with probability one.

If we can only compare two real numbers to within $\epsilon$, estimating $x^\ast$ should be straightforward; you can just use a maximum likelihood method, like you described. If $\epsilon \ll 1$, you should be able to find $x^\ast$ almost exactly from about four samples (on average), and your expected error should be roughly linear in $\epsilon$.

This kind of method should generalize to many probability distributions that are "continuous with a finite number of discrete spikes." The general strategy is to first figure out where the spikes are, using the fact that two equal samples came from a spike with probability 1, and then remove the spikes, leaving you with an ordinary continuous distribution.

-
The zeroes are relevant, it's just that I didn't specify the problem properly (I edited the question). Thanks for the terminology "probability distribution with discrete spikes", that's a good way of describing the problem. The problem I described is a special case of that, in that I know my discrete spike is at 0, and its probability is the cumulative probability of the normal distribution up to -m and the remaining continuous distribution is a Truncated normal distribution. – Emile Apr 14 '10 at 16:49

First, normal random variables are uniquely determined by two moments; e.g., the mean and variance. Here's an easy way if your mean is zero (in which case estimating the parameter $m$ is equivalent to estimating the variance).

Eliminate the zero samples and multiply the non-zero samples by $+1$ with probability $1/2$ and $-1$ otherwise. This will "symmetrize" your sequence; it is an easy exercise to check that the non-zero elements are iid gaussian random variables with mean zero. Then apply your favorite variance estimator to the resulting sequence. If you are confident that they are normal, the best choice is also the unbiased choice, namely $$\hat{ \sigma}^2 = \frac{1}{n}\sum_{i=1}^n g_i^2.$$ Here, $n$ is the number of non-zero elements in your sequence, not the total number of original elements. Also, since you don't have to estimate the mean (because we've assumed that you know it is zero), you need to divide by $1/n$ and not $1/(n-1)$ to get an unbiased estimator.

A final obvious note. You don't actually have to multiply the sequence by $\pm 1$ random variables. You can simply compute the squared sum of your sequence divided by the number of non-zero elements to estimate the variance (equivalent to the equation above). The process I outlined is just a conceptual exercise you can use to justify this computation.

-
Oops, damn, I realized I forgot to add that I generate the data (before truncation) using a distribution with mean m and variance 1. So it's not the variance of the distribution I'm trying to estimate, but the mean. Sorry for the missing information :P – Emile Apr 14 '10 at 16:40

Suppose you perform 2 draws and observe the sequence {v1, 0}. Then your question is:

P( {v1, 0} | m ) = ?.

The intuitive way to think of P( {v1, 0} | m ) is the take the frequentist perspective and imagine that we perform 2 draws N times where N is a very large number and count the number of times we observe {v1, 0}. Then, as N tends to infinity, we have

P( {v1, 0} | m ) = No of times we observe {v1, 0} / N

If you imagine doing this experiment then "No of times we observe {v1, 0} / N" will clearly equal f(v1 | m) * Prob( v < 0 | m ) where v1 ~ N(m,1) I(v1 >0) ) and where v ~ N(m,1).

To see why the above is correct, note that: (a) the two events {v1} and {0} are independent, (b) the number of times you draw v1 is proportional to v1 being drawn from a truncated normal and (c) the number of times you see 0 is proportional to the probability that a normal draw is a negative number.

Estimating m is then a matter of using either maximum-likelihood or bayesian ideas.

Hope that helps.

-

Others have hit on this, but I thought I'd contribute how I'd write the problem down (briefly):

Let $x_i \sim N(m, 1)$ for $i \in \lbrace 1, \dots, N\rbrace$ and define $y_i \equiv x_i$ if $x_i > 0$ and $y_i \equiv 0$ if $x_i \leq 0$. The data in this case is the sample $Y = \lbrace y_1, \dots, y_N \rbrace$. So any likelihood based method (Bayes, MLE, etc) starts with the likelihood for $m$ (or sampling model for $Y$ depending on how you want to look at it). Letting $n$ denote the number of exactly zero observations, this likelihood can be simply expressed as $$\Phi(-m)^n \times \prod_{j \;\mid \; y_j \;\neq \;0} N^+(y_j \;;\; m, 1),$$ where $\Phi(\cdot)$ is the standard Normal CDF and $N^+(\ \cdot \mid \;mean, variance)$ denotes a positive truncated Normal pdf. By explicitly writing the normalizing constant, we can rewrite this as $$\frac{\Phi(-m)^n }{\left(1-\Phi(-m)\right)^{N-n}} \times \prod_{j \;\mid \; y_j \;\neq \;0} N(y_j \;;\; m, 1).$$ This expression makes clear why the zeros $do$ matter, because they are informative about $m$ via the ratio of zeros to non-zeros observed.

-

For the estimate of $m$, I think you can use the nonparametric approach (Nadaraya-Watson estimator) -- it is better than the Bayesian approach. See http://en.wikipedia.org/wiki/Kernel_regression.

-