Added discussion of mean, variance, variance of mean to documentation.

This commit is contained in:
Paul Romano 2012-08-02 12:46:46 -04:00
parent 3c6e800b06
commit 90e093de00

View file

@ -6,29 +6,150 @@ Statistics
As was discussed briefly in :ref:`methods_introduction`, any given result from a
Monte Carlo calculation, colloquially known as a "tally", represents an estimate
of the mean of some random variable of interest. This random variable typically
corresponds to some physical quantity like a reaction rate, a net current across
some surface, or the neutron flux in a region. Given that all tallies are
produced by a stochastic process, there is an associated uncertainty with each
value reported. It is important to understand how the uncertainty is calculated
and what it tells us about our results. To that end, we will introduce a number
of theorems and results from statistics that should shed some light on the
interpretation of uncertainties.
of the mean of some `random variable`_ of interest. This random variable
typically corresponds to some physical quantity like a reaction rate, a net
current across some surface, or the neutron flux in a region. Given that all
tallies are produced by a `stochastic process`_, there is an associated
uncertainty with each value reported. It is important to understand how the
uncertainty is calculated and what it tells us about our results. To that end,
we will introduce a number of theorems and results from statistics that should
shed some light on the interpretation of uncertainties.
--------------------
Law of Large Numbers
--------------------
Let :math:`X_1, X_2, \dots, X_n` be an infinite sequence of independent,
identically-distributed random variables with expected values :math:`E(X_1) =
E(X_2) = \mu`. The sample mean :math:`\bar{X_n} = \frac{X_1 + \dots + X_n}{n}`
converges in probability to the true mean, i.e. for all :math:`\epsilon > 0`
The `law of large numbers`_ is an important statistical result that tells us
that the average value of the result a large number of repeated experiments
should be close to the `expected value`_. Let :math:`X_1, X_2, \dots, X_n` be an
infinite sequence of `independent, identically-distributed random variables`_
with expected values :math:`E(X_1) = E(X_2) = \mu`. One form of the law of large
numbers states that the sample mean :math:`\bar{X_n} = \frac{X_1 + \dots +
X_n}{n}` `converges in probability`_ to the true mean, i.e. for all
:math:`\epsilon > 0`
.. math::
\lim\limits_{n\rightarrow\infty} P \left ( \left | \bar{X}_n - \mu \right |
\ge \epsilon \right ) = 0.
---------------------
Central Limit Theorem
---------------------
------------------------------------------
Estimating Statistics of a Random Variable
------------------------------------------
Mean
----
Given independent samples drawn from a random variable, the sample mean is
simply an estimate of the average value of the random variable. In a Monte Carlo
simulation, the random variable represents physical quantities that we want
tallied. If :math:`X` is the random variable with :math:`N` observations
:math:`x_1, x_2, \dots, x_N`, then an unbiased estimator for the population mean
is the sample mean, defined as
.. math::
:label: sample-mean
\bar{x} = \frac{1}{N} \sum_{i=1}^N x_i.
Variance
--------
The variance of a population indicates how spread out different members of the
population are. For a Monte Carlo simulation, the variance of a tally is a
measure of how precisely we know the tally value, with a lower variance
indicating a higher precision. There are a few different estimators for the
population variance. One of these is the second central moment of the
distribution also known as the biased sample variance:
.. math::
:label: biased-variance
s_N^2 = \frac{1}{N} \sum_{i=1}^N \left ( x_i - \bar{x} \right )^2 = \left (
\frac{1}{N} \sum_{i=1}^N x_i^2 \right ) - \bar{x}^2.
This estimator is biased because its expected value is actually not equal to the
population variance:
.. math::
:label: biased-variance-expectation
E[s_N^2] = \frac{N - 1}{N} \sigma^2
where :math:`\sigma^2` is the actual population variance. As a result, this
estimator should not be used in practice. Instead, one can use `Bessel's
correction`_ to come up with an unbiased sample variance estimator:
.. math::
:label: unbiased-variance
s^2 = \frac{1}{N - 1} \sum_{i=1}^N \left ( x_i - \bar{x} \right )^2 =
\frac{1}{N - 1} \left ( \sum_{i=1}^N x_i^2 - N\bar{x}^2 \right ).
This is the estimator normally used to calculate sample variance. The final form
in equation :eq:`unbiased-variance` is especially suitable for computation since
we do not need to store the values at every realization of the random variable
as the simulation proceeds. Instead, we can simply keep a running sum and sum of
squares of the values at each realization of the random variable and use that to
calulate the variance.
Variance of the Mean
--------------------
The previous sections discussed how to estimate the mean and variance of a
random variable using statistics on a finite sample. However, we are generally
not interested in the *variance of the random variable* itself; we are more
interested in the *variance of the estimated mean*. The sample mean is the
result of our simulation, and the variance of the sample mean will tell us how
confident we should be in our answers.
Fortunately, it is quite easy to estimate the variance of the mean if we are
able to estimate the variance of the random variable. We start with the
observation that if we have a series of uncorrelated random variables, we can
write the variance of their sum as the sum of their variances:
.. math::
:label: bienayme-formula
\text{Var} \left ( \sum_{i=1}^N X_i \right ) = \sum_{i=1}^N \text{Var} \left
( X_i \right )
This result is known as the Bienaymé formula. We can use this result to
determine a formula for the variance of the sample mean. Assuming that the
realizations of our random variable are again identical,
independently-distributed samples, then we have that
.. math::
:label: sample-variance-mean
\text{Var} \left ( \bar{X} \right ) = \text{Var} \left ( \frac{1}{N}
\sum_{i=1}^N X_i \right ) = \frac{1}{N^2} \sum_{i=1}^N \text{Var} \left (
X_i \right ) = \frac{1}{N^2} \left ( N\sigma^2 \right ) =
\frac{\sigma^2}{N}.
We can combine this result with equation :eq:`unbiased-variance` to come up with
an unbiased estimator for the variance of the sample mean:
.. math::
:label: sample-variance-mean-formula
s_{\bar{X}}^2 = \frac{1}{N - 1} \left ( \frac{1}{N} \sum_{i=1}^N x_i^2 -
\bar{x}^2 \right ).
At this point, an important distinction should be made between the estimator for
the variance of the population and the estimator for the variance of the
mean. As the number of realizations increases, the estimated variance of the
population based on equation :eq:`unbiased-variance` will tend to the true
population variance. On the other hand, the estimated variance of the mean will
tend to zero as the number of realizations increases. A practical interpretation
of this is that the longer you run a simulation, the better you know your
results. Therefore, by running a simulation long enough, it is possible to
reduce the stochastic uncertainty to arbitrarily low levels.
------------------------
Random Number Generation
------------------------
@ -85,3 +206,17 @@ so the idea is to determine the new multiplicative and additive constants in
.. _linear congruential generator: http://en.wikipedia.org/wiki/Linear_congruential_generator
.. _Brown: https://laws.lanl.gov/vhosts/mcnp.lanl.gov/pdf_files/anl_rn_arb-strides_1994.pdf
.. _Bessel's correction: http://en.wikipedia.org/wiki/Bessel's_correction
.. _random variable: http://en.wikipedia.org/wiki/Random_variable
.. _stochastic process: http://en.wikipedia.org/wiki/Stochastic_process
.. _independent, identically-distributed random variables: http://en.wikipedia.org/wiki/Independent_and_identically_distributed_random_variables
.. _law of large numbers: http://en.wikipedia.org/wiki/Law_of_large_numbers
.. _expected value: http://en.wikipedia.org/wiki/Expected_value
.. _converges in probability: http://en.wikipedia.org/wiki/Convergence_of_random_variables#Convergence_in_probability