gtMath

We can see that $N_3$ has the lowest variance of the three for all sample sizes. In the legend, I have marked $N_3$ as "MVUE", which stands for minimum variance unbiased estimator. This means that $N_3$ actually has the lowest variance of any unbiased estimator for $N$, not just the three above.

In Part 3 of this post, I will introduce the mathematical heavy machinery (sorry for the pun) to prove that $N_3$ is indeed the MVUE.

Thanks for reading.

Sources

Wikipedia: German Tank Problem
Wikipedia: Discrete uniform distribution
Estimating the Size of a Population (Roger Johnson) - a paper introducing the unbiased estimators above.
An Empirical Approach to Economic Intelligence in World War II (Richard Ruggles and Henry Brodie, 1947) - the original paper describing the use of these statistical methods during WWII.

Since $aI \leq X$ is always true, the means must satisfy the inequality ${\Bbb E}(aI) \leq {\Bbb E}(X) \ ( \spadesuit )$. But the left side of this inequality is equal to: 

$$
\begin{align}
{\Bbb E}(aI) &= a {\Bbb E}(I) \\
&= a(1\cdot {\Bbb P}(I=1) + 0 \cdot {\Bbb P}(I=0)) \\
&= a {\Bbb P}(I=1) \\
&= a {\Bbb P}(X \geq a) \tag{$\dagger$}
\end{align}
$$ Substituting $( \dagger )$ into $( \spadesuit )$ gives
$$ a {\Bbb P}(X \geq a) \leq {\Bbb E}(X) $$ Dividing both sides by $a$ yields the desired inequality.

$\square$

Chebyshev's Inequality: Let $X$ be a random variable with mean $\mu$ and non-zero, finite variance $\sigma^2$, and let $k \geq 1$. Then the probability is at most $\frac{1}{k^2}$ that $X$ falls at least $k$ standard deviations from its mean. In symbols: $$
{\Bbb P}( | X - \mu | \geq k \sigma ) \leq \frac{1}{k^2}
$$ Chebyshev's inequality puts a maximum on how frequently a random variable will exhibit values far away from its mean.

Proof: $|X - \mu| \geq k \sigma$ if and only if $(X-\mu)^2 \geq (k \sigma)^2$. Since $(X-\mu)^2$ is non-negative, we can apply Markov's inequality with $a = (k \sigma)^2$: $$
{\Bbb P}\left[ (X- \mu)^2 \geq (k \sigma)^2 \right] \leq \frac{{\Bbb E}\left[ (X-\mu)^2 \right]}{(k \sigma)^2}
$$ But ${\Bbb E}\left[ (X-\mu)^2 \right]$ is the definition of the variance $\sigma^2$, so $$
{\Bbb P} \left[ (X- \mu)^2 \geq (k \sigma)^2 \right] \leq \frac{\sigma^2}{(k \sigma)^2} = \frac{1}{k^2}
$$ as desired.
$\square$

Now, we can apply Chebyshev's inequality to the sample mean $\overline{X}_n$ to obtain the Weak Law of Large Numbers. For the below, let $X_1, X_2, X_3, \dotsc$ be iid random variables with mean $\mu$ and finite variance $\sigma^2$.

Proof of the (Weak) Law of Large Numbers: In Monte Carlo Simulation - Part 2, we saw that the expected value of the sample mean is the long-run mean, i.e. ${\Bbb E}(\overline{X}_n) = \mu$. Furthermore, the variance of $\overline{X}_n$ is $$
\begin{align}
{\rm var} \left( \overline{X}_n \right) &= {\rm var} \left(\frac{1}{n} \sum_{i=1}^{n}{X_i} \right) \\[2mm]

&= \frac{1}{n^2} \sum_{i=1}^{n}{{\rm var}(X_i)} \\[2mm]

&= \frac{n \sigma^2}{n^2} = \frac{\sigma^2}{n}
\end{align}
$$ Now, let $\epsilon > 0$. Then $$
\begin{align}
k \sigma_{\overline{X}_n} = \epsilon
&\iff k^2 {\rm var} \left( \overline{X}_n \right) = \epsilon^2 \\[2mm]

&\iff \frac{1}{k^2} = \frac{{\rm var}\left( \overline{X}_n \right)}{\epsilon^2} = \frac{\sigma^2}{n \epsilon^2}
\end{align}
$$ Thus, by Chebyshev's inequality, $$
{\Bbb P}\left( \left| \overline{X}_n - \mu \right| \geq \epsilon \right)
\leq \frac{1}{k^2}
= \frac{\sigma^2}{n \epsilon^2}
$$ Taking the limit of both sides yields $$
\lim_{n \rightarrow \infty}{\Bbb P}\left( \left| \overline{X}_n - \mu \right| \geq \epsilon \right)
\leq \lim_{n \rightarrow \infty}{\frac{\sigma^2}{n \epsilon^2}} = 0
$$ Finally, since probabilities are non-negative, the limit must be equal to zero.
$\square$

I hope this post helped clarify what the Law of Large Numbers says and why it's true. Thanks for reading, and feel free to ask any questions in the Comments section.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77

#include<iostream> #include<fstream> #include<random> #include<vector> #include<cmath> #include<chrono> using namespace std; class point //represents a point (x,y) in the plane { private: double x; double y; public: point(double _x, double _y): x(_x), y(_y) {} //constructor void assign_vals(double _x, double _y) {x = _x; y = _y;} double len(void){ //distance from the point to the origin return pow(pow(x,2)+pow(y,2),0.5); } }; // obtain a seed from the system clock: unsigned seed = std::chrono::system_clock::now().time_since_epoch().count(); static std::mt19937 gen(seed); double rand_num(){ static std::uniform_real_distribution<double> urd(-1.0,1.0); //uniform random real number between -1 and 1 return urd(gen); //use gen as a generator } int main() { //drop N random points into the square(-1,1)x(-1,1) //check if they are within circle of radius 1 //circle area = pi*(1^2) = pi //square area = 2*2 = 4 //pi estimate = 4 * (# pts in circle / # pts total) const int N = 100000; //number of sims vector<bool> results(N,0); //vector of N 0's to store intermediate results vector<double> pi_approx_vect(N,0); //vector of N 0's to store approximation results double x, y, r; int sum=0; //sum counter dummy variable double pi_approx; point * p = new point(0,0); //pointer to point ofstream myfile ; myfile.open("Program Output.txt", ios::trunc); for(int j=1; j<=results.size(); j++){ x = rand_num(); //random number btw -1.0 and 1.0 y = rand_num(); //''' p->assign_vals(x,y); r = p->len(); if (r < 1){results[j] = 1;} else {results[j] =0;} sum = sum + results[j]; pi_approx_vect[j] = 4 * sum / static_cast<double>(j); if (myfile.is_open()) { myfile << results[j] << endl; } } delete p; pi_approx = 4 * sum / static_cast<double>(N); cout << "N: " << N << endl; cout << "pi approximation: " << pi_approx; myfile.close(); return 0; }

Any questions or constructive criticism on the code is more than welcome in the Comments section- I am fairly new to C++, so some of the above may be unconventional. Note that it is certainly possible to write the same simulation without using classes, but I did it this way to experiment a bit as a learning exercise for myself.

That will do it for this post. Thanks for reading, and thanks to "Anonymous" for the great Reader Request.

These properties imply the other familiar probability results such as ${\Bbb P}(A^{C}) = 1 - {\Bbb P}(A)$ and, if $A \subseteq B$, ${\Bbb P}(A) \leq {\Bbb P}(B)$.

With its 3 components, a probability space $(\Omega, \mathcal{F}, {\Bbb P})$ serves as the "input" for a model of a random phenomenon. In the next section, we will look at the "output".

Random variables

A random variable (r.v.) is a function$^{\dagger}$ which, given the outcome of a random experiment/state of nature/etc. (represented by an element $\omega \in \Omega$), returns a numerical output of interest. The nomenclature is a bit misleading as a random variable is a function, not a variable. Furthermore, it is not random per se, as the randomness is already captured in the input, $\omega$.

$\dagger$: technically, a random variable must be a measurable function, a restriction meant to exclude functions such as $f(x)=1$ if $x \in A$ and $f(x)=0$ if $x \notin A$, where $A$ is a non-measurable set. While this is a nice example to test the rigor of the theory, it is clearly not the type of function that would ever arise in a practical application, since we would need to invoke the axiom of choice to define a set such as $A$. In short, we don't need to worry about the measurability restriction in practical applications.

As a concrete example, let $\Omega$ represent the set of possible outcomes of rolling five 6-sided dice, so a sample element would be $\omega = (1,4,2,2,6)$. In this case, $\Omega$ is finite with $6^5 = 7{,}776$ elements. Define a function $X: \Omega \rightarrow {\Bbb R}$ where $X(\omega)$ is the sum of the 5 rolls in $\omega$. Note that $X$ is a deterministic (i.e. not random) function based on the random outcome of the dice rolls.

To compute the probability that $X=29$ (for example), we add up the probabilities of the $\omega$'s which will make $X$ equal 29. The only way to have $X=29$ is to have one roll result in 5 and the rest result in 6. Thus, $$
\begin{align}
{\Bbb P}(X=29) &= {\Bbb P}\left( \lbrace \omega \in \Omega: X(\omega)=29 \rbrace \right) \\
&= {\Bbb P}(\text{one 5 and four 6's}) \\
&= (5 \cdot (1/6)) \cdot (1/6)^4 \\
&\approx .06\%
\end{align}
$$ Note for technical readers: by definition, ${\Bbb P}$ takes elements of $\mathcal{F}$ (i.e. subsets of $\Omega$) as inputs and returns a number 0 through 1, but we have written "${\Bbb P}(X=29)$" above. In this slight abuse of notation, we have implicitly used $X$ to "push ${\Bbb P}$ forward" to a probability measure $X_{*}{\Bbb P}$ (the pushforward measure) which takes sets of output values of $X$ as inputs (as opposed to taking subsets of $\Omega$, which actually consist of input values for $X$). Since the interpretation is not ambiguous, it is very common to see an expression such as ${\Bbb P}(X=29)$ written here in place of the more correct $X_{*}{\Bbb P}(\lbrace 29 \rbrace)$.

For a more complicated example of a random variable, let $X$ be the function from $\Omega$ to ${\Bbb R}^{+}$ (positive real numbers- thus, we are ignoring finite minimum tick sizes) which returns the price of some stock at a future time $t$. In this case, an element $\omega \in \Omega$ may represent the history of all trades ever in this stock, which in turn rely on too many phenomena to even imagine. In order to analyze the probability that $X$ takes on a certain value or range of values, we would need to make some assumptions. At the end of this post, we'll see how to do this.

Probability distribution of a random variable

Going back to the example of rolling 5 dice and computing their sum, the possible outputs of $X$ are the integers 5 through 30. Let's call this set of integers $R$ (for the "range" of $X$). We can define a function $f: R \rightarrow [0,1]$ by $f(r) = {\Bbb P}(X=r)$, whose graph would look like this:


Note that $f(r)$ must be between 0 and 1 for each $r$ (i.e. the $y$-axis starts at 0 and goes up to at most 1) since probabilities are always between 0 and 1. Such a function $f$ is called the probability mass function (pmf), or just probability distribution, of the random variable $X$. It is sometimes more convenient to use the cumulative mass function (cmf) of $X$, which is defined by $F(r) = \small{\displaystyle \sum_{j \leq r}}{f(j)}$. From this formula, we can see that $F(r)$ represents the probability that $X$ takes on a value less than or equal to $r$; clearly, $F$ is increasing, i.e. $F(r_1) \leq F(r_2)$ when $r_1 \leq r_2$, and $F(30) = 1$. $f$ and $F$ contain the same information, so we can choose the one which is most convenient for a particular use.

Note also that due to the definitions of $f$ and $F$ as probabilities, it must be the case that $\small{\displaystyle \sum_{r \in R}}{f(r)} = 1$, which is equivalent to the statement above that $F(30) = 1$.

Note: some authors use the term "probability distribution" to refer to the cumulative mass function. In this post and others, I will explicitly use the word "cumulative" if applicable and will use "probability distribution" to refer to the non-cumulative version.

The random variable above is a discrete random variable, meaning it can only take on a countable (in this case, finite) number of values. The stock price example above is a continuous random variable, meaning it can take on an uncountable number of values (in plain English: values from a continuum), in this case ${\Bbb R}^{+}$.

Probability distributions for continuous random variables need to be treated differently since there are so many possible values that the probability of any particular value is effectively zero. Instead, we need to think in terms of ranges of possible values. In this case, we define the probability density function (pdf) $f: {\Bbb R}^{+} \rightarrow [0,1]$ as a function for which $$
{\Bbb P}(a \leq X \leq b) = \int_{a}^{b}{f(x)dx}
$$ Instead of adding the probabilities of values for which $X$ falls between $a$ and $b$, we need to integrate, which is the continuous analog of summation. By analogy to distributions of mass in physics, we also use the word "density" for continuously distributed probability to replace "mass" from the discrete case. However, the term probability density function is more general, and so we often use this term (or its abbreviation, pdf) in the discrete case as well.

Similarly, the cumulative density function (cdf) $F$ would be defined by $$
F(x) = {\Bbb P}(X \leq x) = \int_{- \infty}^{x}{f(t)dt}
$$ where, once again, the summation from the discrete formula has been replaced by integration.

In the stock price example, the pdf of $X$ may look something like this:


As in the discrete case, the probability interpretation of $f$ and $F$ dictates that $$
\int_{-\infty}^{\infty}{f(x)dx} = 1
$$ which is equivalent to $\displaystyle \lim_{x \rightarrow \infty}F(x) = 1$.

When the underlying sample space $\Omega$ is too complicated to model, as in the stock price example, we instead assume that the random variable representing the outcome of interest follows a certain probability distribution. In this way, we "abstract away" the sample space formalism and obtain a concrete tool to analyze the probabilities of possible events. The next section introduces important characteristics of probability distributions, and below, I will introduce a few of the most important probability distributions. Finally, I will explain how we can quantitatively measure whether an assumed probability distribution is reasonable or not.

Mean, variance, and other parameters

Probability distributions come in different shapes and sizes, but there are a few measures available which help characterize them in a standardized way.

The mean of a random variable $X$, typically denoted $\mu$ (or $\mu_X$ when numerous r.v.'s are involved), is essentially the weighted average of the values $X$ can take, weighted by probability of occurrence. Since the mean is also known as the expected value, it is also denoted ${\Bbb E}[X]$. If $f$ is the pdf of $X$, then the mean is calculated as $$
{\Bbb E}[X] = \sum_{j = -\infty}^{\infty}{jf(j)}
$$ for discrete (real-valued) r.v.'s and $$
{\Bbb E}[X] = \int_{-\infty}^{\infty}{x f(x) \, dx}
$$ for continuous (real-valued) r.v.'s.

Getting back to the physics analogy from which the mass/density terminology originates, the mean is the center of mass of the pdf. This means that if the height of the pdf chart represented a weight located at $x$ along the number line, the mean would be the balancing point. Finally, note that the term "expected value" can be a bit misleading, as ${\Bbb E}[X]$ may be a value which $X$ can't actually take on. For example, the "expected" value of a roll of a die is 3.5.

The variance of $X$, usually denoted $\sigma^2$, $\sigma_X^2$, or $\text{Var}[X]$, is the (probability-)weighted average squared distance from the mean. Thus, the formula for variance is $$
\text{Var}[X] = \sum_{j=-\infty}^{\infty}{(j-{\Bbb E}[X])^2 f(j)}
$$ for discrete r.v.'s and $$
\text{Var}[X] = \int_{-\infty}^{\infty}{(x-{\Bbb E}[X])^2 f(x) \, dx}
$$ for continuous r.v.'s. Variance is a measure of the "spread" of $X$: higher variance means that we are more likely to observe values far from the mean. In the mass interpretation, the variance is the moment of inertia: higher variance means it is more difficult to spin the pdf around the mean.

Due to the squaring in the above formulas, variance does not have the same units as $X$ but rather those units squared (e.g. if $X$ is in meters, then $\sigma^2_X$ has units of square meters). For this reason, variance can be cumbersome to interpret, so we often use the standard deviation, defined as $\sigma_X = \sqrt{\sigma^2_X}$, which does have the same units as $X$. Standard deviation measures the average distance from the mean. Note that in order to compute this distance, it is necessary to find the average squared distance and then take the square root: if we tried to directly calculate the expected value of $X - {\Bbb E}[X]$, we would always get zero due to the definition of expected value (since ${\Bbb E}[X]$ is the "balancing point")!

While the mean and variance/standard deviation are useful, they don't capture all the information about a pdf $f$. We could construct many very different pdf's with the same mean and variance, like the two in the figure below:


Other popular measures used to characterize a probability distribution include the median, mode, quartiles, and higher-order moments (i.e. average third, fourth, fifth, etc. powers of distance from the mean). In the interest of length, I won't go into all the different quantitative characteristics of probability distributions, but the take-away point is that in order to fully characterize a random phenomenon, we need more than just the mean and variance.

When someone provides us with a mean statistic, we should typically ask for information on the full probability distribution; the variance/standard deviation is a good start, but we should also look at a histogram, which represents sample data from the distribution: is it unimodal (one "hump") or multimodal (many "humps")? Is it symmetric around the mean? Does it have "fat tails" (i.e. are extreme values likely)? This qualitative information is equally as important as the number crunching for proper statistical analysis.

A few important distributions

In this section, I'll show you three important probability distributions which come up frequently in all sorts of applications.

The first is the uniform distribution, a very simple continuous probability distribution which basically represents "a random number between $a$ and $b$". Since probabilities must add up to 1, i.e. the total area under the curve of the pdf must be 1, the pdf of the uniform distribution is $f(x) = \frac{1}{b-a}$ for all values of $x$ between $a$ and $b$ (and zero otherwise)- in other words, a flat horizontal line. The pdf and cdf are pictured below:


Next is the binomial distribution with parameters $n$ and $p$. This discrete distribution represents the number of "successes" in $n$ random, independent experiments, where each experiment has probability of success $p$ and probability of failure $1-p$ (such experiments are called Bernoulli trials). Since the number of ways to place $k$ successes within $n$ experiments is the binomial coefficient $n \choose k$, and the probability of any particular such sequence is $p^k (1-p)^{n-k}$, the pdf is $f(k) = {n \choose k} p^k (1-p)^{n-k}$. It is relatively simple to prove that the mean and variance of a binomial distribution are $np$ and $np(1-p)$. Below are plots of the pmf and cmf for a few different parameter values:


Last but not least is the Normal distribution (also called Gaussian distribution) with mean $\mu$ and variance $\sigma^2$, the pervasive "bell curve", which has pdf's and cdf's as shown below for a few parameter values:


As you can see from the pdf diagram, larger values of $\mu$ move the peak to the right, while larger values of $\sigma$ make the distribution fatter (and thus shorter, as the total area underneath must always be 1). At first glance, the formula for the Normal pdf, $$
f(x) = \frac{1}{\sigma \sqrt{2 \pi}} e^{-\frac{(x-\mu)^2}{2\sigma^2}}
$$ looks daunting. However, thanks to various favorable properties of the exponential function (especially certain integral formulas involving exponentials in the integrand), calculations with the Normal distribution are actually quite tractable. In addition, this distribution conveniently provides error bands based on the 68-95-99.7 rule: ~68% of data from a Normal distribution falls within 1 standard deviation ($1\sigma$) of the mean, ~95% of the data falls within 2 standard deviations ($2\sigma$), and ~99.7% falls within 3 standard deviations ($3\sigma$):


Thus, if we reasonably believe that a certain random variable has a Normal distribution, and we calculate a mean $\bar{x}$ and standard deviation $s$ from sample data, then we can assume that data will fall within $2s$ of $\bar{x}$ 95% of the time, a very simple yet powerful tool for confidence interval analysis.

The Normal distribution applies well to certain natural phenomena such as people's heights, errors in measurements, or the $x$-coordinate of the position of a particle diffusing in a liquid$^{\ddagger}$. However, in the next post on Monte Carlo methods, we will see that the bell curve also arises as the distribution of the sample means when we repeat a random experiment many times, a result known as the Central Limit Theorem.

$\ddagger$: The fact that this is Normal is part of the mathematical definition of Brownian motion. Towards the end of the linked post, you will see (1) how the Normal distribution arises in the context of diffusion and (2) how Gaussian integral formulas make Normal distribution math a piece of cake.

There are many other well-known probability distributions, but hopefully the three above give you a sense of what they can look like. In the final section of this post, we'll look briefly into how to choose the right distribution to model a random outcome of interest.

Testing for goodness of fit

Above, I mentioned that we can ignore impossibly complex sample spaces by simply assuming a particular probability distribution for a random variable we'd like to study. The question is, which probability distribution(s) can we reasonably assume? While a thorough treatment is beyond the scope of this post, at a high level, the process, known as hypothesis testing, typically proceeds as follows:

1. Hypothesize, i.e. guess, a distribution for the random variable in question. This is called the null hypothesis. Often, the null hypothesis is that the data follows a Normal distribution.
2. Gather sample data and calibrate the distribution using the sample mean, standard deviation, etc.
3. Assuming the random variable follows the (calibrated) hypothesized distribution, calculate how likely it is that hypothetical data from a sample of a certain size would differ from the theoretical distribution by at least as much as your sample data did. This is called a goodness-of-fit test. There are numerous well-known statistical goodness-of-fit tests such as the $\chi^2$ ("chi-squared") test.
4. The probability calculated in step 3 is called the $P$-value, and if it is lower than a pre-set threshold (usually 5% or 1%), i.e. if it's very unlikely that the sample data came from the hypothesized distribution, then we conclude that we must reject the null hypothesis.

In this post, we saw how the technical underpinnings of probability theory make their way into the definition of random variables, which model outcomes of random phenomena. We saw that a random variable is actually a function, not a variable, and that its input, not its output, drives the randomness in the model. Since the inputs may be far too complicated to describe directly, we can instead assume a certain probability distribution for the output and test our assumption using a statistical goodness-of-fit test. We often test for Normality so that we can take advantage of the 68-95-99.7 rule or other favorable properties to simplify our analysis.

As a simple example, suppose you're an investment banker, and you want to lend a client money overnight, taking 100,000 shares of a stock from his portfolio as collateral. If the client fails to repay, you will liquidate the stock tomorrow at the then-prevailing price. If the stock trades at USD 100.00 per share right now, and you lend the client USD 9,800,000.00, is the 2% buffer enough to protect you in the unlikely case the client fails to repay?

You download 5 years of historical data and plot the following histogram of the stock's daily % returns:

The mean of the sample data is +0.06% and the standard deviation is 0.92%, so you decide to test the null hypothesis that this stock's daily returns are Normally distributed with mean 0.06% and standard deviation 0.92%. You frantically search gtmath.com for how to actually do a $\chi^2$ test and find that it's not on there; luckily, it's all over Google, and you can even do it in Excel. You run the test and obtain a $P$-value of 0.42. Since this is not less than 5%, you conclude that the Normality assumption was not unreasonable, and thus even if the client fails to repay, the probability is still less than 2.5% (half of 5% since we only care about downward moves) that the stock drops by more than 1.84% by tomorrow, so you deem the 2% buffer sufficient and go ahead with the loan.

In the next post, I will demystify the buzzword "Monte Carlo simulation" and also introduce the Law of Large Numbers and Central Limit Theorem. For now, thanks for reading!