General Definitions

The Probability Mass Function (PMF) is another name for the Probability Distribution Function (PDF), $P(x)$, of a discrete random variable X expresses the probability that $X$ takes the value $x$. I.e, $$P(x)\equiv P(X=x),\mathrm{\forall }x$$ The PMF has the following properties: $$\begin{gathered} 0\le P(x)\le 1\mathrm{\forall }x \\ \sum _{x}P(x)=1 \end{gathered}$$

The Cumulative Mass Function (CMF) or Cumulative Probability Function (CPF), $F(x_0)$, for a random variable $X$, gives the probability that X does not exceed the value $x_0$: $$F(x_0)=P(X\le x_0)=\sum _{x\le x_0}P(x)$$ The CMF has the following properties: $$\begin{gathered} 0\le F(X_0)\le 1,\mathrm{\forall }{x}_0 \\ {x}_0<x_1⟹F(x_0)<F(x_1) \end{gathered}$$

Expected Value and Variance

Expected value defined as: $$E(X)=\mu =\sum _{x}xP(x)$$ Where $\mu$ is called the mean and is the mean of the population.

The calculation of measurements like $E(X^2)$ becomes... $$E(X^2)=\sum _x{x}^{2}P(x)$$

Variance, ${\sigma }^2$ is the average squared distance from the mean and is defined as follows. The square is taken so that distances don't cancel eachother out (i.e, a negative distance and positive distance could result is a very small average distance, which is not what we want). $$\begin{array}{rl}\text{Var}(X)={\sigma }^2& =\sum _x(x-\mu {)}^{2}P(x)\\ & =E[(X-\mu {)}^2]\end{array}$$ We can re-write this as follows: $$\begin{array}{rl}\text{Var}(X)={\sigma }^2& =\sum _x(x-\mu {)}^{2}P(x)\\ & =\sum _x(x^2-2x\mu +{\mu }^2)P(x)\\ & =\sum _x{x}^{2}P(x)-2\mu \sum xP(x)+{\mu }^2\sum P(x)\\ & =E[X^2]-2\mu E[X]+E[X{]}^2\\ & =E[X^2]-2E[X{]}^2+E[X{]}^2\text{ , (because }\mu ==E[X]\text{ and }\therefore \mu E[X]==E[X{]}^2\text{)}\\ & =E[X^2]-E[X{]}^2\end{array}$$ In summary, $$\text{Var}(X)=E[(X-\mu {)}^2]=E(X^2)-{\mu }^2=E(X^2)-E(X{)}^2$$

Standard Deviation, $\sigma$, is the positive square root of the variance.

The calculation of measurements like $\text{Var}(aX)$ becomes... $$\begin{array}{rl}\text{Var}(aX)& =E[(aX{)}^2]-E[aX{]}^2\\ & =E[a^2{X}^2]-(aE[X]{)}^2\\ & =a^{2}E[X^2]-a^{2}E[X{]}^2\\ & =a^2(E[X^2]-E[X{]}^2)\\ & =a^2\text{Var}(X)\end{array}$$

Why does $E[(X-\mu {)}^2]$ equal $E(X^2)-{\mu }^2$? $$\begin{array}{rl}E[(X-\mu {)}^2]& ={\mathrm{\Sigma }}_{i}P()[x_i-\mu {]}^2\\ & ={\mathrm{\Sigma }}_{i}P(x_i)[x_i^2-2x_i\mu +{\mu }^2]\\ & ={\mathrm{\Sigma }}_{i}P(x_i)x_i^2-2\mu {\mathrm{\Sigma }}_{i}P(x_i)x_i+{\mu }^2{\mathrm{\Sigma }}_{i}P(x_i)\end{array}$$ From our definition of $E(X)$ we know that ... $${\mathrm{\Sigma }}_{i}P(x_i)x_i^2=E(X^2)$$ And... $$\begin{array}{rl}2\mu {\mathrm{\Sigma }}_{i}P(x_i)x_i& =2\mu E(X)\\ & =2E(X)E(X)\\ & =2E(X{)}^2\end{array}$$ We also know that the sum of all the probabilities in the distribution must sum to 1, so therefore... $${\mu }^2{\mathrm{\Sigma }}_{i}P(x_i)={\mu }^2$$ Thus we can say this... $$\begin{array}{rl}E[(X-\mu {)}^2]& =E(X^2)-2E(X{)}^2+E(X{)}^2\\ & =E(X^2)-E(X{)}^2\end{array}$$ Yay!

For example, lets say that we have weighted dice so that the probabilities are as follows:

The population mean, $\mu$ becomes: $$E[X]=\mu =\sum _{x}xP(x)=1\times 0.05+0.1\times 2+0.1\times 3+0.1\times 4+0.1\times 5+0.55\times 6=4.75$$ Of course, the dice does not have a face value of 4.75, but over many many rolls this would be the average score. The variance and standard deviation are calculated in the same way using the above formulas.

Linear Functions Of X

The expected value and variance of a linear function of $X$ is another fairly useful result and we will use it to get some very important results in the section on sampling distributions. For example, we could arbitrarily define a function over dice rolls. Not sure why we'd do this, but say we said the experiement result was $g(X)=2X+5$ where $X$ is the random variable which gives the dice face rolled. Now we have:

The expected value of $g(x)$ is therefore: $$E[g(X)]=\mu =\sum _{x}xP(x)=7\times 0.05+0.1\times 9+0.1\times 11+0.1\times 13+0.1\times 15+0.55\times 17=14.5$$ Interestingly $14.5=2\times 4.75+5$! It looks like $E[g(X)]=g({\mu }_x)$, and for linear functions this is the case (but not so if $g(X)$ is not linear!). When $g(X)$ is LINEAR and $g(X)=a+bX$: $$\begin{array}{rl}{\mu }_{g(X)}& =E[a+bX]\\ & =a+b{\mu }_X\\ \\ {\sigma }_{g(X)}^2& =var[a+bX]\\ & =b^2{\sigma }_X^2\end{array}$$

We can derive the first equation as follows, using our first definition of the expected value of a random variable. $$\begin{array}{rl}E[a+bX]& =\sum (a+b{x}_i)P(x_i)\\ & =\sum aP(x_i)+\sum b{x}_{i}P(x_i)\\ & =a+b\sum x_{i}P(x_i)\\ & =a+bE[X]\\ & =a+b{\mu }_X\end{array}$$ The same can be done for the variance can also be done in a similar manner. TODO: add this here.

Binomial Distribution

Bernoulli Model

Experiment with two mutually exclusive and exhaustive outcomes. One has the probability $p$ and the other has probability $(1-p)$. $$\begin{array}{rl}P(0)& =P(fail)\\ & =(1-p)\\ \\ P(1)& =P(success)\\ & =p\end{array}$$ Therefore, using the formulas for mean and variance from the previous sections, we can say the following. $$\begin{array}{rl}{\mu }_x& =E[X]\\ & =\sum _{x}xP(x)\\ & =(0)(1-p)+(1)p\\ & =p\\ \\ {\sigma }_x^2& =E[(X-\mu {)}^2]\\ & =\sum (x-{\mu }_x{)}^{2}P(x)\\ & =p(1-p)\end{array}$$

Binomial Distribution

Bernoulli experiment repeated $n$ times where the repetitions are independent. This is like selection with replacement.

We know from previous discussions that if two events are independent then $P(A\cap B)=P(A)P(B)$ and by extension that $P(A\cap B\cap \cdots \cap Z)=P(A)P(B)\dots P(Z)$. Therefore if I do $n$ experiments and want to know the probablity of $m$ successes in a row and then $n-m$ failures in a row, the probability is: $$\begin{array}{rl}P(msuccesses)& =P(1)P(1)\dots P(1)P(0)P(0)\dots P(0)\\ & =p^m(1-p{)}^{n-m}\end{array}$$ Where there are $m$ $P(1)$'s in the above equation and $(n-m)$ $P(0)$'s. BUT this would be the probability of getting $m$ successes in a row and then $(n-m)$ failures in a row. The question, however, doesn't care about the specific order: we dont care if we get SSFF.. or SFSF... or SSFS... and so on. We need to figure out how many of these combinations there are and account for this!

Let's take a simple example. I have a bag with 5 balls in it: 3 blue, 2 red. What is the probability of drawing 1 blue ball and 1 red ball if I select using replacement (to make the selections independent). Replacement means that what ever ball I pick first, I record the result and then put it back in the bag before making my next pick. In this way, the probability of selecting a particular colour does NOT change per pick.

Note that there are only 2 colours of ball in our bag... this is because we are talking about an experiement where there are only two outcomes, labeled "success" and "failure". We could view selecting a blue as "success" and a red as "failure", or vice versa.

So... selecting a blue and a red ball. Sounds like $P(\text{blue}\cap \text{red})$ right? Well, almost, but not quite. Take a look at the sample space below, out of which, the events of interested are highlighted in a light blue colour.

Here I am using an alphabetical subscript to indicate the specific ball. I.e. $B_x$ is a different ball to $B_y$. The order of selection is given by the order of writing. I.e., "$B_x{B}_y$" means $B_x$ was picked on the first turn and then $B_y$ was picked on the second turn.

There are two clear groups of outcomes that will satisfy the question. We see that in one group we drew a red ball first and in the other we drew a blue ball first. So, there are $5\times 5$ total possible events, and of these $6+6$ are of interest. Therefore, $$\begin{array}{rl}P(\text{drawing a blue and a red})& =\frac{6}{25}+\frac{6}{25}\\ & =2\times \frac{6}{25}\end{array}$$ But, hang on minute! Isn't $P(\text{drawing a blue and a red})$ the same as $P(B\cap R)$?! Well, no, as we can see below... $$P(B\cap R)=P(B)\times P(R)=\frac{3}{5}\times \frac{2}{5}=\frac{6}{25}$$ So the two expressions are clearly not the same thing. To clear up my confusion I asked the guys at Maths Exchange, and got the following awesome answer from a very nice chap named Joriki. Im quoting it (almost) verbatim because it was just so good!

Note, that in Jorki's example $R_1$ and $R_2$ do not refer to different balls: the subscripts refer to different draws. Therefore, $R_1$ and $R_2$ could be the same red ball drawn on turn $1$ and turn $2$.

So, having understood this, we can see that to get the total probability of $m$ successess (1 blue ball) out of $n$ trials (2 selections), in the case where events are independent, we are concerned with the number of combintations in which the events can occur.

The formula for the binomial distribution is... $$\begin{array}{rl}P(m\text{successes in }n\text{ independent trials})& =P(m)\\ & ={}_n{C}_m\cdot p^m\cdot (1-p{)}^{(n-m)}\\ & =\frac{n!}{m!(n-m)!}\cdot p^m\cdot (1-p{)}^{(n-m)}\end{array}$$ For $m=0,1,\dots ,n$.

This, incidentally is the same as asking for $P((n-m)\text{failures in }n\text{ independent trials})$ because $$\begin{array}{rl}{}_n{C}_{(n-m)}& =\frac{n!}{(n-(n-m))!(n-m)!}\\ & =\frac{n!}{m!(n-m)!}\\ & ={}_n{C}_m\end{array}$$

Let's do this to exhaustion... let's imagine another bag. It doesn't matter how many blue and red balls are in there, just so long as I can select 3 blue balls and 1 red. $P(\text{blue})=p$ and $P(\text{red})=(1-p)$. Once we've made this selection of 3 blues, 1 red, the question is, how many ways were there of getting to this outcome. Now let's check our understanding of using combinations (esp. vs. permutations) to get this...

We can draw a little outcome tree as follows:

Clearly there are 4 ways to arrive at the selection of 3 blues and 1 red, once we have made that particular selection. Remember we're not asking for a blue/red on a specific turn... we just care about the final outcome, irrespective of the order in which the balls were picked. And we can see... $$\begin{array}{rlr}{}_4{C}_3& =\frac{4!}{3!(4-3)!}& =4\end{array}$$

Why don't we use permutations? The answer is that we don't care if we got $B_a{B}_b{B}_c{R}_a$ or $B_b{B}_a{B}_c{R}_a$ or $B_c{R}_a{B}_a{B}_b$ etc etc, where, as before, the alphabetical substrcipts distinguish the ball, not the turn on which it is drawn as that is given by order of writing.

Having gone through this we can then use the formulas for expected value, or mean, and variance. To make the notation similar to previous examples where we used $P(x)$, in the above formula I just change $m$ for $x$ so that we get $$P(x)=\frac{n!}{x!(n-x)!}\cdot p^x\cdot (1-p{)}^{(n-x)}$$ Where $x=0,1,\dots ,n$, meanding that this is the probability of $x$ successes out of $n$ trials.

Recalling that $$E(X)=\mu =\sum _{x}xP(x)$$ We can say that the population mean for the Bernoulli distribution is $$\begin{array}{rl}E(X)=\mu & =\sum _{x=0}^{x\le n}x\cdot (\frac{n!}{x!(n-x)!}\cdot p^x\cdot (1-p{)}^{(n-x)})\end{array}$$ And... the rest of the proof gets a little complicated... this PDF by Joel Feldman gives the derivation, which I found by doing a quick google. Liked the explanation.

The population mean and variance are summarised as follows. $$\begin{array}{rl}\mu & =E[X]\\ & =nP\\ \\ {\sigma }_x^2& =E[(X-\mu {)}^2]\\ & =nP(1-P)\end{array}$$

Often when talking about a binomial distribution you will see something like $$P(X|n=?,p=?)$$ This is the binomial distribution for $X$ successes out of $n$ trials with the probability of success given by $p$.

We can plot some example distributions... lets do this using Python.

import matplotlib.pyplot as pl
import numpy as np
from matplotlib.font_manager import FontProperties
from scipy.stats import binom

fontP = FontProperties()
fontP.set_size('small')

n=50
pLst=[0.1, 0.25, 0.5, 0.75, 0.9]
x = np.arange(-1, n+2)
fig, ax = pl.subplots()

for p in pLst:
   dist = binom(n, p)
   ax.plot(x, dist.pmf(x),linestyle='steps-mid')

ax.legend(['$p={}$'.format(p) for p in pLst],
          ncol = 3,
          prop=fontP,
          bbox_to_anchor=[0.5, 1.0],
          loc='upper center')
fig.show()
fig.savefig('binomial_distrib.png', format='png')

Poisson Distribution

Didn't like the starting explanation in [1] so had a look on Wikipedia and found a link to the UMass Amherst Uni's stats page on the Poisson distribution which I though was really well written. That is the main reference here.

Poisson distribution gives the probability of an event occuring some number of times in a specific interval. This could be time, distance, whatever. What matters is that the interval is specific and fixed.

The example used on the UMass Amherst is letters received in a day. The interval here is one day. The poisson distribution will then tell you the probability of getting a certain number of letters in one day, the interval. Other examples could include the number of planes landing at an airport in a day or the number of linux server crashes in a year etc...

The interval is one component. The other is an already observed average rate per interval, or expected number of successes in an interval, $\lambda$. For example, we might have observed we get on average 5 letters per day ($\lambda =5$), or that 1003 planes land at our airport ($\lambda =1003$) per day or that there are 4 linux server crashes per year ($\lambda =4$) etc...

So, poisson has an interval and an observed average count for that interval. The following assumptions are made:

• The #occurrences can be counted as an integer
• The average #occurrences is known
• The probability of the occurence of an event is constant for all subintervals. E.g., if we divided the day into minutes, the probability of receiving a letter in any minute of the day is the same as for any other minute.
• There can be no more than one occurrence in the subinterval
• Occurrences are independent

The distribution is defined as follows, where $\lambda$ is the expected number of events per interval. $$P(X=x)=\frac{e^{-\lambda }{\lambda }^x}{x!},x\ge 0$$

Eek... didn't like the look of trying to derive the mean and variance. The population mean and variance are as follows. $$\begin{array}{rl}{\mu }_x& =E(X)=\lambda \\ \\ {\sigma }_x^2& =E[(X-{\mu }_x{)}^2]=\lambda \end{array}$$

We can plot some example distributions... lets do this using Python.

import matplotlib.pyplot as pl
import numpy as np
from matplotlib.font_manager import FontProperties
from scipy.stats import poisson

fontP = FontProperties()
fontP.set_size('small')

expectedNumberOfSuccessesLst = [1, 5, 10, 15]
x = np.arange(-1, 31)
fig, ax = pl.subplots()

for numSuccesses in expectedNumberOfSuccessesLst:
   ax.plot(x, poisson.pmf(x, numSuccesses),linestyle='steps-mid')

ax.legend(['$\lambda={}$'.format(n) for n in expectedNumberOfSuccessesLst],
          ncol = 3,
          prop=fontP,
          bbox_to_anchor=[0.5, 1.0],
          loc='upper center')
fig.show()
pl.show()
fig.savefig('poisson_distrib.png', format='png')

Lets take a little example. Looking at the UK government report on road casualties in 2014, there were a reported 194,477 casualties in 2014. This gives us an average of 532.8137 (4dp) casualties per day! So, we could ask some questions...

What is the probability that no accidents occur on any given day? Recalling that... $$P(X=x)=\frac{e^{-\lambda }{\lambda }^x}{x!}$$ We will set $\lambda =\frac{194,477}{365}$ and $x=0$ and plug these into the above forumla... $$\begin{array}{rl}P(X=0|\lambda =\frac{194,477}{365})& =\frac{e^{(-194,477/365)}{\left(\frac{194,477}{365}\right)}^0}{0!}\\ & =\frac{e^{(-194,477/365)}\times 1}{1}\\ & =3.998994e-232\end{array}$$ Okay... so the probability that on any given day in the UK that there are no road casualties is pretty (very scarily) small!

Of course, we could critisize this analysis as not taking into account seaons, or weather conditions etc etc. In winter, for example, when the roads are icy one might expect the probability of an accident to be greater, but for the purposes of a little example, I've kept it simple.

We might also ask, what is the probability that the number of casualties is less than, say, 300? $$P(X<300)=\sum _{x=0}^{299}P(X=x)$$ Eek! Don't really want to be calculating this by hand so let's use Python's SciPi package:

from scipy.stats import poisson
print poisson.cdf(299, 194477.0/365)

This outputs $2.24102710469e-28$ which is still pretty unlikely!. For $X<400$ we get $9.60741455596e-10$, for $X<500$ we get $0.0783182362497$. So we can imagine that the distribution is really quite steep:

import matplotlib.pyplot as pl
import numpy as np
from scipy.stats import poisson

x = np.arange(200, 800)
fig, ax = pl.subplots()
ax.plot(x, poisson.pmf(x, 194477/365), linestyle='steps-mid')
ax.set_title('Probability of #road casualties per day in UK')
fig.show()
fig.savefig('poisson_road_fatilities_distrib.png', format='png')

Joint Distributions

When one thing is likely to vary in dependence on another thing. For example, the likely longevity of my milk will correspond to some extent to the temperature of my fridge: there is a relationship between these two variables so it is important that any model produced includes the effect of this relationship. We need to define our probabilities that random variables simultaneously take some values.

Enter the joint probability function. It expresses the probability that each random variable takes on a value as a function of those variables.

A two variable joint probability distribution is defined as: $$P(x,y)=P(X=x\cap Y=y)$$ This is a subtle difference in the notation of $P(x,y)$ from our previous examples in the starting sections where we wrote things like $P(A\cap B)$. The former is a distribution where the random variables simultaneously take on the values $x$ and $y$ respectively and the latter is just the probability of two individual events happening simultaneously.

The marginal probabilities are the probabilities that one random variable takes on a value regardless of what the other(s) are doing. In our two variable distribution we have: $$\begin{array}{rl}P(x)& =\sum _{y}P(x,y)\\ \\ P(y)& =\sum _{x}P(x,y)\end{array}$$

Joint probability functions have these properties:

1. $0\le P(x,y)\le 1$ for any pair of values $x$ and $y$,
2. The sum of the join probabilities is 1.

The conditional probability function looks like this... $$\begin{array}{r}P(y|x)=\frac{P(x,y)}{P(x)}\\ \\ P(x|y)=\frac{P(x,y)}{P(y)}\end{array}$$

The random variables are independent if their joint probability function is the product of their marginal probability functions: $$P(x,y)=P(x)P(y)⟹\text{independence}$$ When they are independent we also have $$\begin{array}{rl}\text{independence}& ⟹P(y|x)=P(y)\text{and}\\ & ⟹P(x|y)=P(x)\end{array}$$

The conditional mean and variance are: $$\begin{gathered} {\mu }_{Y|X}=E[Y|X]=\sum _{Y}yP(y|x)\text{and} \\ {\mu }_{X|Y}=E[X|Y]=\sum _{X}xP(x|y) \end{gathered}$$

Covariance and Correllation

Measure of joint variablility: the nature and strength of a relationship between two variables.

Covariance between two random variables $X$ and $Y$ is $Cov(X,Y)$, given by: $$\begin{array}{rl}Cov(X,Y)& =E[(X-{\mu }_x)(Y-{\mu }_y)]\\ & =\sum _x\sum _y(x-{\mu }_x)(y-{\mu }_y)P(x,y)\end{array}$$ This is equivalent to the following: $$\begin{array}{rl}Cov(X,Y)& =E[XY]-{\mu }_x-{\mu }_y\\ & =\sum _x\sum _{y}xyP(x,y)-{\mu }_x{\mu }_y\end{array}$$ A covariance that is strongly negative indicates a good inverse linear relationship. Strongly positive indicates a good linear relationship and near zero indicates no relationship. Note that if two random variables are statistically independent, the covariance between them is (near) zero. But, if the covariance is (near) zero it does not necessarily mean there is no relationship, it might just not be linear.

Correlation is just a "normalisation" of the covariance such that the measure is limited to being in the range [-1, 1]. $$Corr(X,Y)=\frac{Cov(X,Y)}{{\sigma }_x{\sigma }_y}$$

This YouTube video explains the difference really well:

Cumulative Distribution Function (CDF)

CDF gives probability that continuous random variable $X$ does not exceed a particular value, $x$: $$CDF(x)=P(X<x)$$ The probability that $X$ takes on a signular value is zero. I.e., $P(X==x)=0$. Therefore, it doesn't matter if we write $<$ or $\le$...

To get the probability that $X$ lies within a range use: $$\begin{array}{rl}P(x_1<X<x_2)& =CDF(x_2)-CDF(x_1)\\ & =P(X<x_2)-P(X<x_1)\end{array}$$

Probability Density Function (PDF)

Call our PDF $f(x)$... Then the total area under the curve $f(x)$ is 1. I.e., $${\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}f(x)\mathrm{d}x=1.0$$ The probability that $X$ takes a value between two limits is... $$P(x_1<X<x_2)={\int }_{x_1}^{x_2}f(x)\mathrm{d}x$$ And the CDF is also defined by.... $$CDF(x_0)={\int }_{min(X)}^{x_0}f(x)\mathrm{d}x$$ We can judge whether a distribution is a valid PDF using the above definitions.

Uniform Distribution

The uniform distribution is one where the probability that $X$ takes a value in a fixed size interval is the same regardless of where the interval "starts" (as long as the interval is contained entirely within the range for which $X$ is non-zero). I.e, $$f(x)=\{\begin{array}{ll}\frac{1}{x_u-x_l},& x_l\le X\le x_u\\ 0,& \text{otherwise}\end{array}$$

Normal Distribution

This distribution is defined by a gaussian bell-shaped curve. The probability distribution function for a normally distributed random vairable, $X$, is given by this equation: $$f(X=x)=\frac{1}{\sqrt{2\pi {\sigma }^2}}\cdot e^{-(x-\mu {)}^2/2{\sigma }^2}$$ For any normal distribution the following applies:

1. About 68% of the data will fall within one standard deviation, $\sigma$, of the mean, $\mu$,
2. About 95% of the data will fall within two standard deviations, $2\sigma$, of the mean, $\mu$,
3. Over 99% of the data will fall within three standard deviations, $3\sigma$, of the mean, $\mu$,
4. Mean is $E[X]=\mu$,
5. Variance is $Var(X)=E[(X-\mu {)}^2]={\sigma }^2$,

The normal distribution is defined with the following notation: $$X\sim N(\mu ,{\sigma }^2)$$ There is no simple algebraic expression for the cummulative distribution function. Can't really say I fancy the idea of integrating the above function! There are many numerical approximations. Computer could do it for us, but another way is to convert every normal distribution to the standard normal distribution.

The standard normal distribution is a normal distribution where the the mean is 0 and the variance is 1: $$Z\sim N(0,1)$$ And $$Z=\frac{X-\mu }{\sigma }$$ In English we can say that the z-score/value is the number of standard deviations $X$ is away from the (population) mean $\mu$.

To solve $P(a<X<b)$ we solve the following: $$P(\frac{a-\mu }{\sigma }<Z<\frac{b-\mu }{\sigma })$$ For which we can consult a standard normal table.

We can also get the probability of a Z score and vice versa in Python by using scipy.stats.norm.cdf() and scipy.stats.norm.ppf() respectively as documented here.

Percentiles

Wikipedia defines a percentile as a measure giving the value below which a given percentage of obersvations fall. If an observation is AT the x^th percentile then it has the value, V, at which x% of the scores can be found. If it is IN this percentile then its value is in the range 0 to V.

Percentiles aren't really specific to normal distributions but you often get asked questions like "what is the 95^th percentile of this distribution?"

In the standard normal distribution...

1. 90^th percentile is at z=1.28. (This is because $P(z<1.28)=90$).
2. 95^th percentile is at z=1.645. (This is because $P(z<1.645)=95$).
3. And so on...

A Little Intro...

So far we have been talking about population statistics. The values $\mu$ and $\sigma$ have been the mean and standard deviation of the population. However, generally it is pretty impossible to gather information about an entire population: this can be due to the cost that would be involved, or perhaps that time that such an endevour would take, for example. It might also be undesirable to analyse an entire population if, for example, analysis involved destruction of the samples taken!

So, what is normally done is to take a sample from the population and then use the sample statistics to make inferences about the population statistics. The image below shows three samples that have been taken from a population. Each sample set can, and will most probably, have a different shape, mean, and variance!

We can demonstrate this concept using a quick little Python program to take 4 samples from the normal distribution where each sample has 10 members:

import numpy as np
import matplotlib.pyplot as pl

numSampleSets = 4
numSamplesPerSet = 10

# Limit individual plots, otherwise the first plot takes forever and is rubbish
doShowIndividualSamples = (numSampleSets <= 8) and (numSamplesPerSet <= 50)

# Create a numSampleSets (rows) x numSamplesPerSet (cols) array where we use
# each row as a sample.
randomSamples = np.random.randn(numSampleSets, numSamplesPerSet)

# Take the mean of the rows, i.e. the mean of each sample set
means = randomSamples.mean(axis=1)

if doShowIndividualSamples:
   fig, axs = pl.subplots(nrows=int(numSampleSets/2.0+0.5), ncols=2)
   xticks = np.arange(numSamplesPerSet, dtype='float64')
   for idx in range(numSampleSets):
         ax_col = idx % 2
         ax_row = int(idx/2.0)
         thisAx = axs[ax_row][ax_col]
         thisAx.bar(xticks, randomSamples[idx][:], width=1)
         thisAx.set_xticks(xticks + 0.5)
         thisAx.set_xticklabels(xticks, fontsize=8)
         thisAx.axhline(y=means[idx], color="red", linewidth=2)
         thisAx.set_title("Sample set #{}".format(idx))
         thisAx.grid()

xticks = np.arange(numSampleSets, dtype='float64')
fig2, ax2 = pl.subplots(nrows=2)
ax2[0].bar(xticks, means)
ax2[0].grid()
ax2[0].set_title("Distribution of {} sample means (n={})".format(numSampleSets, numSamplesPerSet))
ax2[0].set_ylabel("Mean value")
ax2[0].set_xlabel("Sample set #")

ax2[1].hist(means, 50)
ax2[1].grid()
ax2[1].set_title("Histogram of {} sample means (n={})".format(numSampleSets, numSamplesPerSet))
ax2[1].set_ylabel("# sample mean's")
ax2[1].set_xlabel("Sample mean bin")

pl.tight_layout()
pl.show()

The script above produces the following graphs. The x-axis is just the sequence number of the sample member and the y-axis the value of the sample member. The horizontal line is the sample mean.

We can see from this little example that the samples in each of the 4 instances are different and that the mean is different for each sample. Keep in mind that the graph shown when you run the above script will be different as it is a random sample :)

So we can see that although the population mean is fixed the sample mean can vary from sample to sample. Therefore, the mean of a sample is itself a random variable and as a random variable it will have it's own distribution (same applies for variance).

But, we want to use the sample statistics to infer the population statistics. How can we do this if the sample mean (and variance) can vary from sample to sample? There are a few key statistical theories that will help us out...

Independent and Identically Distributed (I.I.D.) Random Variables

Wikipedia says that ...a sequence or other collection of random variables is independent and identically distributed (I.I.D.) if each random variable has the same probability distribution as the others and all are mutually independent...

It is often the case that if the population is large enough and a sample from the population only represents a "small" fraction of the population that a set of simple random samples without replacement still qualifies I.I.D. selections. If you sample from a population and the sample size represents a significant fraction of the population then you cannot assume this to be true.

We're going to need to know this a little later on when we talk about conistentent estimators and the law of large numbers...

The Law of Large Numbers

The law of large numbers states that the average of a very large number (note we haven't quite defined what "very large" is) of items will tend towards the average of the popultion from which the items are drawn. I.e., the sample mean (or any other parameter) tends towards the population mean (or parameter, in general) as the number of items in the sample tends to infinity.

This is shown empirically in the example below:

import numpy as np
import matplotlib.pyplot as pl

maxSampleSizes = [100, 1000, 100000]
fig, axs = pl.subplots(nrows = 3)
for idx, ssize in enumerate(maxSampleSizes):
   sample_sizes = np.arange(1,ssize+1)
   sample_means = np.random.randn(ssize).cumsum() / sample_sizes
   axs[idx].plot(sample_sizes, sample_means)
   axs[idx].set_xlabel("Sample size, max is {}".format(ssize));
   axs[idx].set_ylabel("Average")
   axs[idx].axhline(y=0, color="red")

fig.show(); pl.show()

Explain the code a little... np.random.randn(ssize) generates a numpy array with ssize elements. Each element is a "draw" from the standard normal distribution. The function cumsum() then produces an array of ssize elements where the second element is the sum of the first 2 elements, the third element is the sum of the first 3 elements and so on. Thus we have an array where the i^th element is the sum of i samples. Dividing this by the array sample_sizes gives us an array of means where the i^th element is the mean of a sample with i items.

Running the code produced the figure below...

To summarise, the law of large numbers states that the sample mean of I.I.D. samples is consistent with the population mean. The same is true for sample variance. I.e., as the number of items in the sample increases indefinitely, the sample mean will become a better and better estimate of the population mean.

Sampling Distribution Of The Sample Mean

We've established that between samples, the mean and variance of the samples, well... varies! This implies that we can make a distribution of sample means.

The "sampling distribution of the sample mean" (a bit of a mouthful!) is the probability distribution of the sample means obtained from all possible samples of the same number of observations.

That's a bit of a mounthful! What it means is that if we took, from the population, the exhaustive set of all distinct samples of size $n$, and then took the mean of each sample, we could figure out what the propability is that any sample of size $n$ has a specific mean. Thus we build up the probability distribution of sample means, where the sample size is $n$. We will see that this probability distribution is centered around the mean of all the sample means and that it also has a normal distribution (see LLN and CLT).

For example if I am a factory owner and I produce machines that accurately measure distance, I could have a population of millions. I clearly do not want to test each and every device coming off my production line, especially if the time that testing requires is anything approaching the time taken to produce the item: I'd be halving (or more) my production throughput!

So what do I do? I can take a sample of say 50 devices each day. I can test these to see if they accurately measure distance, and if they do I can assume the production process is still running correctly.

But as we have seem if I test 50 devices each day for 10 days, each of my 10 sample sets will have a different mean accuracy. On the first day, the mean accuracy of 50 devices might be 95%, on the second day, the mean of the next 50 devices might be 96.53% and so on.

The sample mean has a distribution. As we can take many samples from a population we have sampled the sample mean, so to speak, hence the rather verbose title "sampling distribution of the sample mean".

Now, I can ask, on any given day, "what is the probability that the mean accuracy of my 50 devices is, say 95%?". This is the sampling distribution of the sample mean: the probability distribution of the sample means (mean accuracy of a sample of devices) obtained from all possible samples (theoretical: we can't actually measure all possible samples!) of the same number of observations (50 in this case).

The law of large numbers gives us a little clue as to how variable the sample means will be... we know that if the sample size is large then the sample mean tends towards the population mean as has a much narrower variance about the population mean. So, 50 samples is better than 10. But 10,000 samples would be even better. The question is how many samples do we need. Will try to answer that (much) later on...

Sample Mean And It's Expected Value (The Mean Of Sample Means)

Continuing with the example of distance measurement devices. We could say each device in our sample set is represented by random variables $A_1,A_2,...,A_{50}$ (we samples 50 devices each day). We can generalise this to $m$ devices in the sample set. Using our definition of expected value we can define the sample mean as the random variable ${\bar{A}}_i$, for the i^th sample set of, and is defined as follows. $${\bar{X}}_i=\frac{1}{m}\sum _{i=1}^m{A}_i$$

Now we can calculate the expected value of the sample mean, $E(\bar{X})$ (the mean of the sample means).

Imagine that we take $n$ samples, where each sample is a set of $m$ devices from our population. If we represent the mean of each sample as the random variables ${\bar{X}}_1,{\bar{X}}_2,\cdots ,{\bar{X}}_n$, then the expected value of the distribution of sample means is the mean expected value of the samples... $$E(\bar{X})=E(\frac{1}{n}({\bar{X}}_1+{\bar{X}}_2+\cdots +{\bar{X}}_n))$$ We saw that expected value of a linear combination of random variables is the linear combination of each random variable's expected value. This means that we can write the following. $$E(\bar{X})=\frac{1}{n}[E({\bar{X}}_1)+E({\bar{X}}_2)+\cdots +E({\bar{X}}_n)]$$ The law of large numbers also tells us that if the size of each sample is sufficiently large, it's sample mean, $\overline{X_i}$ will tend towards the population mean, $\mu$. So, if $m$ is large enough we can say... $$\begin{array}{rl}E(\bar{X})& =\frac{1}{n}[\mu +\mu +\cdots +\mu ]\\ & =\frac{n\mu }{n}\\ & =\mu \end{array}$$for the This means that the mean of "the sample distribution of sample means" is the population mean when the size of each sample ($m$) is sufficiently large. Put another way, the distribution of samples means in centred around the population mean.

A single sample mean can therefore be larger or smaller than the population mean, but on average, there isn't any reason to expect that it is either. As the sample size increases we also know that the likelihood of the sample mean being higher or lower than the population mean decreases (law of large numbers).

The Variance Of The Sample Mean w.r.t. Sample Size And Standard Error

Sample variance can be written as ${\sigma }_{\bar{X}}^2$ or $S^2$. $$\text{Var}(\bar{X})={\sigma }_{\bar{X}}^2=S^2=\frac{{\sigma }^2}{n}$$ We can see that variance of the sampling distribution decreases as sample size increases. Therefore the larger the sample the more accurately we can infer population statistics.

How did we get the above results? Well, it all is a bit brain melting, but here goes. The variance of the mean of sample means is written as $\text{Var}(\bar{X})$. Because it is the variance of the mean of sample means, we can write this: $$\text{Var}(\bar{X})=\text{Var}\left(\frac{1}{n}\sum _{i=1}^{i\le n}{\bar{X}}_i\right)$$ The constant $n$ can be taken out of the variance. TODO - finish this off.

The sample standard deviation is often refered to as the standard error. $$\text{Standard error}={\sigma }_{\bar{X}}=\frac{\sigma }{\sqrt{n}}$$ If sample size not super small compared to population use finite population correction factor $\frac{N-n}{N-1}$ to get $$\text{Var}(\bar{X})=\frac{{\sigma }^2}{n}\cdot \frac{N-n}{N-1}$$

The Shape Of The Samplign Distribution Of Sample Means

When samples are selected from a population that has a normal distribution we will see (anecdotally) that the sampling distribution of sample means is also normally distributed.

A histogram of the sample means gives a better picture of the distribution of sample means from the 4 sample sets. Diagrams have been drawn using the previous python example.

Now observe what happens when we increase the number of sample sets taken to 50. The distribution of sample means is still pretty distributed without much definition of the shape of the distribution...

Now observe what happens when we increase the number of sample sets taken to 1000. The distribution of sample means begins to take shape...

Now observe what happens when we increase the number of sample sets taken to 1000. The distribution of sample means now really looks very much like a gaussian distribution... interesting!

Clearly when we sample from the population the sample can be, to varying degrees, either a good or bad representation of the population. To help ensure that that sample represents the population (i.e., no "section" of the population is over or under represented) a simple random sample is usually taken.

Also of great interest in the above is that as the number of samples taken increases the distribution of sample means appears to become more and more gaussian in shape!. What is always true, and even more important is that the sampling distribution becomes concentrated closer to the population mean as the sample size increases, as we have seen in the maths earlier. The example above also doesn't proove it... its just anecdotal evidence. We could also notice that as #samples increases the variance of the sample mean decreases, again as we saw in the maths earlier.

Remember: take care to differentiate between the sample mean, $\bar{x}$, sometimes witten as ${\mu }_{\bar{X}}$, and the population mean, $\mu$. We do not know the population mean. We know the sample mean. The population mean is fixed. The sample mean, as seen, will vary between samples.

We must also remember to differentiate between the sample variance, ${\sigma }_{\bar{X}}^2$ and the population variance, ${\sigma }^2$. We do not know the population variance. We know the sample mean. The population variance is fixed. The sample variance, as seen, will vary between samples.

Also remember that the samples were drawn from a population that is normally distributed. So the sampling distribution was also normally distributed, as we have anacdotally seen. The central limit theory (CLT) will show us that is doesn't matter what the population distriution is. Enough samples and the sampling distribution of the sample mean will always tend towards normal distribution.

Because the sampling distribution of sample means is normally distributed, it means that we can also form a standard normal distribution for the sample means, which we do as follows: $$\begin{array}{rl}Z& =\frac{\bar{X}-\mu }{{\sigma }_{\bar{X}}}\\ & =\frac{\bar{X}-\mu }{\frac{\sigma }{\sqrt{n}}}\end{array}$$