Average from distribution-fit vs sample average

Question

If I know the functional form of a distribution, is it better to fit the distribution to the data to find the mean, or is it better to just look at sample average.

I know sample average is the unbiased estimator of the mean, but I feel like fitting to a known distribution may reduce the effect of outliers. Is that true?

If not, is there a way to use the information about the shape of distribution to find a better estimation of the mean? What about other moments?

Edit By fitting the distribution to the data, what I had in mind was something like FindDistributionParameter function in Mathematica. This is what it does:

The maximum likelihood method attempts to maximize the log-likelihood function $\sum\log(f(x_i;\theta))$, where $\theta$ are the distribution parameters and $f(x_i;\theta)$ is the PDF of the symbolic distribution.

Answer 1

There's a few issues here.

First, for a few popular distributions (eg Normal, Binomial, Poisson, Gamma), the maximum likelihood estimator of the mean is just the sample average, so it won't make any difference

For symmetric distributions with sensible parametrisations, the maximum likelihood estimator of the mean and the sample average will both estimate the same quantity. If the distribution you're assuming when you do the estimation is correct, the maximum likelihood estimator will be more precise. The advantage can be quite large for heavy-tailed distributions. If the distribution you're assuming is not correct (but everything is still symmetric), then it's not obvious which estimator will be better.

A good example here is the double-exponential distribution. The MLE of the mean is the sample median, which is more precise than the sample average if the data really come from this distribution or something close, but if the data really come from a Normal distribution the sample median is less precise than the sample average. A large fraction of the field of robust statistics studied this sort of question.

Once you stop assuming symmetry, it's still true that the MLE is more accurate than the sample average when the assumed distribution is correct, but when the assumed distribution is incorrect the two estimators aren't even estimating the same thing. Take the double-exponential distribution again. If you assume that's the distribution (so that the MLE is the sample median), but the data really come from a Gamma distribution or some other asymmetric distribution then the sample median no longer estimates the same thing as the sample mean. It's a lot harder to say anything about robust estimation, because different estimators typically estimate different things.

The sample average has the good feature that it always estimates the mean (if one exists). No other estimator has this feature. The cost is that if you know the true distribution well enough you can get a more accurate estimate than the sample average.

And finally, because it's easy to confuse the sample average and the distribution mean, a lot of people try not to use the term "mean" for other estimators of the mean. If you had a double-exponential distribution and you worked out the maximum likelihood estimator of the mean, you'd probably want to call it the median so that people didn't get confused and think you meant the sample average.

Answer 2

Take a look at the MATLAB code below:

clear; close all; rng(1)
ntrials = 1000;
nsamples = 10000;

%% beta distribution: true mean 0.9

xbetamu = nan(ntrials, 1);
xbetapmu = nan(ntrials, 1);
for i = 1:ntrials
    xbeta = betarnd(90, 10, nsamples, 1);
    xbetap = betapdf(xbeta, 90, 10);

    xbetamu(i) = mean(xbeta);
    xbetapmu(i) = sum(xbeta.*xbetap)/sum(xbetap);
end

figure; hist([xbetamu xbetapmu], 100)
title('beta distribution')
legend('sample mean', 'pdf mean')

%% gamma distribution: true mean 4

xgammamu = nan(ntrials, 1);
xgammapmu = nan(ntrials, 1);
for i = 1:ntrials
    xgamma = gamrnd(2, 2, nsamples, 1);
    xgammap = gampdf(xgamma, 2, 2);

    xgammamu(i) = mean(xgamma);
    xgammapmu(i) = sum(xgamma.*xgammap)/sum(xgammap);
end

figure; hist([xgammamu xgammapmu], 100)
title('gamma distribution')
legend('sample mean', 'pdf mean')

%% normal: true mean 5

xnormmu = nan(ntrials, 1);
xnormpmu = nan(ntrials, 1);
for i = 1:ntrials
    xnorm = normrnd(5, 2, nsamples, 1);
    xnormp = normpdf(xnorm, 5, 2);

    xnormmu(i) = mean(xnorm);
    xnormpmu(i) = sum(xnorm.*xnormp)/sum(xnormp);
end

figure; hist([xnormmu xnormpmu], 100)
title('normal distribution')
legend('sample mean', 'pdf mean')

%% poisson: true mean 5

xpoismu = nan(ntrials, 1);
xpoispmu = nan(ntrials, 1);
for i = 1:ntrials
    xpois = poissrnd(5, nsamples, 1);
    xpoisp = poisspdf(xpois, 5);

    xpoismu(i) = mean(xpois);
    xpoispmu(i) = sum(xpois.*xpoisp)/sum(xpoisp);
end

figure; hist([xpoismu xpoispmu], 100)
title('poisson distribution')
legend('sample mean', 'pdf mean')

For any of the distributions with nonzero skewness, the distribution of the pdf-weighted means is shifted in the direction of the skewness with respect to the distribution of the sample mean. It also looks like it has a smaller variance. Here, we used the true parameter values of the distributions; it seems like fitting the distributions would only perform as well or worse since the fitting process itself requires using the samples.

The statement that the sample mean is an unbiased estimator of the mean holds regardless of the distribution of the data: \begin{align} \mathbb{E}[\bar{x}] = \frac{1}{n}\sum_{i=1}^n\mathbb{E}[x_i] = \frac{1}{n}\sum_{i=1}^n\theta = \frac{n\theta}{n} = \theta \end{align} Since this doesn't depend on the distribution of the data, your estimate will always be unbiased.

On the other hand, if you have outliers, there's a chance you could grossly mis-estimate your mean if you don't have knowledge of the functional form of the distribution (see code and figures below).

%% normal with outliers: true mean somewhere above 5

xnormmu = nan(ntrials, 1);
xnormpmu = nan(ntrials, 1);
for i = 1:ntrials
    xnorm = normrnd(5, 2, nsamples, 1);
    if logical(binornd(1, 0.1, 1))
        xnorm = [xnorm; 20000];
    else
        xnorm = [xnorm; normrnd(5, 2, 1, 1)];
    end    

    xnormp = normpdf(xnorm, 5, 2);

    xnormmu(i) = mean(xnorm);
    xnormpmu(i) = sum(xnorm.*xnormp)/sum(xnormp);
end

figure; hist([xnormmu xnormpmu], 100)
title('normal distribution with outliers')
legend('sample mean', 'pdf mean')

%% poisson with outliers: true mean somewhere above 5

xpoismu = nan(ntrials, 1);
xpoispmu = nan(ntrials, 1);
for i = 1:ntrials
    xpois = poissrnd(5, nsamples, 1);
    if logical(binornd(1, 0.1, 1))
        xpois = [xpois; 20000];
    else
        xpois = [xpois; poissrnd(5, 1, 1)];
    end

    xpoisp = [poisspdf(xpois, 5)];

    xpoismu(i) = mean(xpois);
    xpoispmu(i) = sum(xpois.*xpoisp)/sum(xpoisp);
end

figure; hist([xpoismu xpoispmu], 100)
title('poisson distribution with outliers')
legend('sample mean', 'pdf mean')

Here, 0.1 is the probability of having an outlier in your dataset. Approximately 10% of the sample-mean estimates of the population mean are much farther from the true mean than the pdf-weighted estimates. In this case, knowledge of the true distribution more frequently results in a better estimate of the population mean than just using the sample mean. The sample mean is still ``unbiased'' (because you have to weight the means of the datasets with and without outliers to get the true mean, which is no longer just 5 in either of the two above cases). However, since you know that outliers exist and that their values are inflating the mean (which is actually 5), the sample estimate is much worse on average compared to the pdf-fitted mean.