Normal distribution construction in Python

Here is an example with the normal distribution that will seem trivial after the t-distribution (here).

The basic form of the normal is exp {-x2/2}. We define that as a Python function f(x), vectorize it, and construct an array X of discrete points from -10 to +10 with interval dx = 0.001. We apply the vectorized function to the array to get the relative densities. So that we obtain the correct area under the curve, we multiply the height (value of f(x)) of each piece by its width, dx. When we sum up all the pieces, the total is equal to √2π as seen in the printout. We divide by this value to normalize the distribution so that its total area is equal to 1 and it becomes a pdf.

The form of the normal that includes a term for the standard deviation is normal(x). Everything is as before, except we substitute z for x, and at the end we find that the normalizing constant is 1/σ√2π. We plot it to have something pretty to look at.

[UPDATE:
In my discussions of probability distributions of late I've played a little fast and loose with terminology. The pdf (probability density function) has a value for any x. At x = μ the exponential term is equal to 1, and the value is 1/σ√2π. For the standard normal, this is about 0.4.

>>> 1.0/sqrt(2*pi) 0.3989422804014327

Of course, this is somewhat misleading, since the probability that x = any particular point approaches zero, because x is a real number with infinitely many values. The way we actually use the pdf is to ask what is the probability for a particular window, a range of values for a < x < b, and at least conceptually, this is done by integrating the function between these limits. That's where the discrete version that I developed comes in handy. We evaluate p(x) for a small enough interval and then multiply by the interval size to get a probability for that slice. To integrate between limits, just add the included slices. But for this to work, we have to "normalize" the pdf so that the total of all the little rectangles of width dx is equal to 1. It also helps to have the discrete version in developing the cdf (cumulative distribution function), since we can just accumulate the pdf as we move along the values of X generating the cdf as we go.

Technically, the pdf(x) is the slope of the cdf(x), which gets around the issue mentioned above for a continuous function. ]

output:

2.507 2.507 5.013 5.013

code listing:

from __future__ import division import numpy as np import math import matplotlib.pyplot as plt @np.vectorize def f(x): return math.e**(-0.5*(x**2)) @np.vectorize def normal(x,mu=0,sigma=1): z = (x-mu)/sigma return math.e**(-0.5*(z**2)) #================================== dx = 0.001 X = np.arange(-10,10+dx,dx) pdf = f(X) pdf *= dx print round(sum(pdf),3), S = math.sqrt(2*math.pi) print round(S,3), pdf /= S sigma = 2 pdfn = normal(X,0,sigma) pdfn *= dx print round(sum(pdfn),3), Sn = math.sqrt(2*math.pi) * sigma print round(Sn,3), pdfn /= Sn #================================== plt.plot(X,pdf,color='r',lw=4) plt.plot(X,pdfn,color='k',lw=4) ax = plt.axes() ax.set_xlim(-6,6) m = max(pdf) ymin,ymax = -m/100,m*1.1 ax.set_ylim(ymin,ymax) plt.text(3,0.75*ymax,s='$\sigma = 1$', color='r',fontsize=24) plt.text(3,0.65*ymax,s='$\sigma = 2$', color='k',fontsize=24) plt.savefig('example.png')

More on sampling a distribution

Some time ago I had a post about sampling from a distribution. More specifically, the problem is to generate random samples given a cdf (cumulative distribution function). There have been a couple of useful comments, and I'd like to extend the explanation as well.

Previously, a function was defined that can be used to generate a probability for any value of a random variable x, from a probability density function or distribution (pdf) with a given mean and standard deviation. We remember that for a continuous distribution, the actual probability at any discrete x is zero, since the total number of possible x's is infinite. Technically the definition is that the pdf is the derivative of the cdf, which I misremembered as the cumulative density function. So the cdf(x) is the area under the pdf from negative infinity to x.

One nice thing about this approach is it is then easy to define a sum of weighted distributions.

We get what looks like a smooth curve by plotting a (relatively) large number of points (in this example, 3502). The plotting uses matplotlib (see this post referencing set-up on the Mac). The cdf is computed by simply accumulating values from the pdf. Normalization is usually done by dividing by the total, but the method I showed was just slightly more subtle:

cdf *= dx

.

Unfortunately, it is also wrong! (Sorry). I won't try to explain what I was thinking, but the fact that it gave an accurate value (as shown by the maximum of the cdf being equal to 1) is an accident, and you should instead simply divide by the sum of the values, and moreover, do the operation on the pdf before constructing the cdf:

pdf /= sum(pdf)

This change to the pdf means we need to magnify it before plotting, and really, should provide a different y-axis on the right hand side, with the true values. The left-hand y-axis only is accurate only for the cdf.

The idea for sampling is to generate a random float using np.random.random(), and then ask which of the values in the cdf (which by definition are ordered), first exceeds this value. The indexes resulting from repeating this procedure are concentrated in the steep part of the cdf (the peaks of the pdf), because the probability that a given position in our "discretized" form of the cdf satisfies this relationship is proportional to the slope of the cdf curve (the added vertical distance between a given index and the one previous).

As a reader suggested, an improvement to the code is to recognize that the list we're searching (the cdf) is ordered and so can more efficiently be searched using a binary search. The code is a little tricky to write, so I skipped it last time. And luckily, Python comes with "batteries included" and for this application what we want is the bisect module from the standard library. The example find_le function's docstring says: 'Find rightmost value less than or equal to x'. I just modified this code (which calls bisect.bisect_right) to return the index rather than the value.

We're interested in intervals where the slope is steepest. The original find_first function returned the index of the right-hand value, while this function will return the index of the left-hand value. I suppose either one is fine, but perhaps it would be better to use the midpoint.

There are a few more steps in the code that are a little obscure, including bins with fractional width, and the use of the Counter class to organize the data for the histogram. But this post is getting a bit long so I'll skip them for now.



Modified code:

import math, sys, bisect import numpy as np import matplotlib.pyplot as plt import Counter def normal(mu,sigma): def f(x): z = 1.0*(x-mu)/sigma e = math.e**(-0.5*z**2) C = math.sqrt(2*math.pi)*sigma return 1.0*e/C return f p1 = normal(0,2) p2 = normal(10,1) p3 = normal(18,0.5) # sum of weighted normal distributions @np.vectorize def p(x): return 0.5*p1(x) + 0.25*p2(x) + 0.25*p3(x) dx = 0.01 xmax = 25 R = np.arange(-10,xmax+dx,dx) # dashed lines plt.plot((R[0],R[-1]),(0,0),color='k',ls=':',lw=2) plt.plot((R[0],R[-1]),(1,1),color='k',ls=':',lw=2) plt.plot((R[0],R[-1]),(-0.5,-0.5),color='k',lw=2) pdf = p(R) S = sum(pdf) pdf /= S plt.plot(R,pdf*(len(R)/10.0),color='r',lw=3) #print len(R) S = sum(pdf) print S #=========================================== cdf = [pdf[0]] for e in pdf[1:]: cdf.append(cdf[-1] + e) cdf = np.array(cdf) #cdf /= S plt.plot(R,cdf,color='b',lw=3) ax = plt.axes() ax.set_xlim(-6,xmax) ax.set_ylim(-0.55,1.05) #=========================================== def find_le(n,L): 'modified from bisect.find_le' i = bisect.bisect_right(L,n) if i: return i raise ValueError samples = list() width = 0.2 f = 1/width for i in range(10000): n = np.random.random() # must adjust to actual range value = find_le(n,cdf)*dx - 10.0 # trick to truncate at fractional values samples.append(np.floor(f*value)/f) #=========================================== c = Counter.Counter(samples) maxn = c.most_common(1)[0][1] for k in c: n = c[k] n = n * 0.45 / maxn r = plt.Rectangle((k,-0.5), width=width,height=n, facecolor='magenta') ax.add_patch(r) plt.savefig('example.png')

Sampling a distribution

I was listening to an introductory talk about MCMC (link). I'll have more to say about the talk another time (hint: I think it's terrific). The speaker described a method for sampling from a probability distribution. The idea is to generate random samples from U[0,1], and then find the value of x at which the cumulative density function first exceeds the random number. The resulting samples of x are plotted as a histogram.

The Python script has three parts. In the first section we set up a weighted mixture of three normal distributions of different mean and sd, and plot that in red using x increments of 0.01. In the second section the same values are used to generate a discrete cdf for the same points. This is plotted in blue (after normalizing by the interval size).

Then we implement the algorithm for sampling. Binning for the histogram used the collections.Counter class that will be in Python 2.7. I don't have that, so I grabbed the class from here.

You might notice that the labels on the y-axis also cover the histogram. It would be better to have two separate plots, but I'm not quite sure how to handle that yet.

[UPDATE: more here]

import math, sys import numpy as np import matplotlib.pyplot as plt import Counter def normal(mu,sigma): def f(x): z = 1.0*(x-mu)/sigma e = math.e**(-0.5*z**2) C = math.sqrt(2*math.pi)*sigma return 1.0*e/C return f p1 = normal(0,2) p2 = normal(10,1) p3 = normal(18,0.5) # sum of weighted normal distributions @np.vectorize def p(x): return 0.5*p1(x) + 0.25*p2(x) + 0.25*p3(x) dx = 0.01 xmax = 25 R = np.arange(-10,xmax+dx,dx) # dashed lines plt.plot((R[0],R[-1]),(0,0),color='k',ls=':',lw=2) plt.plot((R[0],R[-1]),(1,1),color='k',ls=':',lw=2) plt.plot((R[0],R[-1]),(-0.5,-0.5),color='k',ls=':',lw=2) L = p(R) plt.plot(R,L,color='r',lw=3) #=========================================== cdf = [L[0]] for e in L[1:]: cdf.append(cdf[-1] + e) cdf = np.array(cdf) cdf *= dx plt.plot(R,cdf,color='b',lw=3) ax = plt.axes() ax.set_xlim(-6,xmax) ax.set_ylim(-0.55,1.05) #=========================================== def find_first(n,L): for i,e in enumerate(L): if n < e: return i return len(L) samples = list() width = 0.4 f = 1/width for i in range(1000): n = np.random.random() # must adjust to actual range value = find_first(n,cdf)*dx - 10.0 # trick to truncate at fractional values samples.append(np.floor(f*value)/f) #samples = np.array(samples) c = Counter.Counter(samples) maxn = c.most_common(1)[0][1] for k in c: n = c[k] n = n * 0.45 / maxn r = plt.Rectangle((k,-0.5), width=width,height=n, facecolor='green') ax.add_patch(r) plt.savefig('example.png')

Gamma distribution

I posted previously about the beta distribution here, here and here. We ran into it because it is the conjugate prior for Bayesian analysis of problems in binomial proportion. That's because the likelihood function is in the same form.

The beta distribution is a continuous probability distribution with two shape parameters α and β (or a and b). If we consider p as the random variable (and q = 1-p), then the unnormalized distribution is just:

f(x) = pa-1 qb-1

The function has symmetry: if we switch a and b, and also p and q, everything would look the same. Varying the values for a and b can yield a wide variety of shapes. For large a and b, the plot approaches a normal distribution with mean = a / a+b.

The need to normalize the beta distribution brings us to the gamma distribution. This one can be viewed in a simple way but also can be fairly complex. The simple version is that, for positive integers, Γ(x) is just (x-1)!.

The normalizing constant for the beta distribution with parameters a and b is:

Γ(a+b) / (Γ(a) * Γ(b))

The gamma distribution is frequently used in phylogenetic models, for example, to model the distribution of variability during evolutionary time among different positions in a protein. The gamma distribution also has two parameters but these apparently have quite different roles.

They are called the shape parameter (k) and the scale parameter (θ). More generally, the unnormalized gamma distribution is:

f(x) = xk-1 exp { -x / θ }

We can get an idea of the roles of the two variables by considering this:

mean = kθ variance = kθ2

Let's look at some plots obtained using different values for the two parameters. In the series below, k goes from 1 to 4, and then for each plot theta varies from 0.25, 0.5, 1, 2, 3 as the color goes from blue to green to magenta to red to maroon.

N = 6 thetaL = c(0.25,0.5,1,2,3) color.list = c('blue','darkgreen', 'magenta','red','darkred') #k = 1,2,3,4 plot(1:N,ylim=c(0,0.8),type='n') for (i in 1:5) { curve(dgamma(x, shape=k,scale=thetaL[i]), from=0,to=N,lwd=3, col=color.list[i],add=T) }

Geometric distribution

According to wikipedia, the geometric distribution is the probability distribution which describes the number of Bernoulli trials required to obtain a single success.

So, in the simple case of a fair coin (p = 1/2):

P(X=1) = 1/2 P(X=2) = 1/4 (1/2 times 1/2) P(X=3) = 1/8 (1/2 times 1/4)

According to mathworld, the geometric distribution is the only discrete memoryless random distribution, and is a discrete analog of the exponential distribution.

The memoryless property can be seen easily if we take the geometric series:

1/2 + 1/4 + 1/8 ...

If we have already obtained a failure on the first trial, then we remove the first term (corresponding to success on the first trial) and then normalize by dividing by the sum of all remaining terms (1 - 1/2):

= 2/4 + 2/8 + 2/16 ... = 1/2 + 1/4 + 1/8 ...

The normalization is needed because we require the sum of all the terms to add up to 1 for a proper probability distribution. In general, for process with probability of success p and failure q = 1 - p:

P(X=k) = qk-1 * p

According to wikipedia, the mean is 1/p and the variance is q/p². We can compare that to the exponential distribution with a pdf of:

λe-λx

and mean = 1/λ and variance = 1/λ².

It is not clear to me at present why the expressions for the variance don't match up.

R code:

x=numeric(7) x[1] = 1 for (i in 2:length(x)) { x[i] = x[i-1]/2 } plot(x,type='s', ylim=c(0,1.0), col='red',lwd=3)

Examining distributions in R

In this post, I want to look at some distributions using R. Let's get 100 samples from the normal distribution with mean = 5:

set.seed(157) v = rnorm(100,mean=5) summary(v)

> summary(v) Min. 1st Qu. Median 2.375 4.292 5.032 Mean 3rd Qu. Max. 5.072 5.865 8.349

The sort function does what it says. We make a plot called a "step" plot:

n = length(v) s = sort(v) plot(s,(1:n)/n, type='s',lwd=2,col='darkred', ylim=c(0,1)) points(mean(v),0.5, pch=16,cex=2,col='blue')

You can also use draw a Q-Q plot (quantile-quantile). This plots the quantiles for our distribution (y-axis) against those for a normal distribution (x-axis). A straight line indicates that our distribution is close to normal.

qqnorm(v)

Let's get a larger sample:

v = rnorm(10000,mean=5) m = mean(v)

> m [1] 4.988964

s = sd(v)

> s [1] 1.000791

The quantile function accepts an argument telling it how to make the quantiles:

S = seq(0,1,by=0.001) x = quantile(v,S) round(x[25:27],2) round(x[975:977],2)

> round(x[25:27],2) 2.4% 2.5% 2.6% 3.04 3.05 3.06 > round(x[975:977],2) 97.4% 97.5% 97.6% 6.93 6.95 6.98

As expected, since the mean is 5 and the standard deviation is 1, 97.5 % of the values are < 5 + 1.96.

I showed this trick for plotting the normal distribution before. Now we plot it in red, and then overlay it with the density from the t-distribution (obtained with the function dt) with degrees of freedom (df) = 9,6,3,2,1.

plot(function(x) dnorm(x), -5:5,max(v),lwd=2,col='red') color.list=c( 'purple','blue','steelblue', 'darkred','magenta') L = c(9,6,3,2,1) for (i in 1:5) { x = L[i] plot(function(x) dt(x,df=i), -5:5,max(v),lwd=2, col=color.list[i], add=T) }

The t distribution has fatter tails. If we had plotted it for larger values of df, it would asymptotically approach the normal distribution.

By using a very large sample, we can get a good idea for the cut-offs for 95% and 97.5% from the normal distribution:

S = seq(0,1,by=0.001) v = rnorm(100000) x = quantile(v,S) round(x[950:952],2) round(x[975:977],2)

> round(x[950:952],2) 94.9% 95.0% 95.1% 1.63 1.64 1.65 > round(x[975:977],2) 97.4% 97.5% 97.6% 1.94 1.96 1.98

And now compare them to the cut-offs for the t-distribution with df=9. We'll use these in the next example.

w = rt(100000,df=8) y = quantile(w,S) round(y[950:952],2) round(y[975:977],2)

> w = rt(100000,df=8) > y = quantile(w,S) > round(y[950:952],2) 94.9% 95.0% 95.1% 1.84 1.85 1.86 > round(y[975:977],2) 97.4% 97.5% 97.6% 2.28 2.30 2.33

Student's t test

This post is a simple demo of using R to carry out Student's t test.

Let's look at a population of values with a normal distribution, mean = 5 and standard deviation = 1.

set.seed(157) v = rnorm(10000,5,1)

We draw 4 samples without replacement:

e1 = sample(v,4)

> round(e1,2) [1] 5.34 4.80 5.56 4.96

A t test for one sample tests the null hypothesis that the mean μ for the population from which the sample is drawn is equal to μ₀. For example it could be that we have many observations of untreated cells (from which we get μ₀), and now we wish to estimate whether the mean values of treated cells are detectably different.

result = t.test(e1,mu=6)

The argument alternative = 'two.sided' is the default, so we don't need to specify it here.

> result One Sample t-test data: e1 t = -4.7613, df = 3, p-value = 0.01759 alternative hypothesis: true mean is not equal to 6 95 percent confidence interval: 4.608507 5.723439 sample estimates: mean of x 5.165973

Even with only four samples and a difference in means of (6-5) / 6 the result of the t test tells us that we can reject the null hypothesis that μ = μ₀ = 6, with p=0.018.

Now, it might have been the case that before we saw the data (and that proviso is crucial), we expected from the nature of the treatment that the mean of treated population would be less than the untreated population. In that case, we would be justified in specifying a one-sided test:

result = t.test(e1,mu=6,alternative='less')

We note the p-value is:

p-value = 0.008796

In reality, because of biological variation (as well as unintended variation in experiment conditions) we would always include a control group for such an experiment.

w = rnorm(10000,6,1) e2 = sample(w,4) result = t.test(e1,e2,alternative='less')

> result Welch Two Sample t-test data: e1 and e2 t = -1.3385, df = 3.416, p-value = 0.1314 alternative hypothesis: true difference in means is less than 0 95 percent confidence interval: -Inf 0.6194634 sample estimates: mean of x mean of y 5.165973 6.084764

Now it is much more difficult to see a result with significance. If there 25 samples in the control group we can see a difference:

result = t.test(e1, rnorm(25,6,1),alternative='less')

> result Welch Two Sample t-test data: e1 and rnorm(25, 6, 1) t = -3.4116, df = 11.547, p-value = 0.002717 alternative hypothesis: true difference in means is less than 0 95 percent confidence interval: -Inf -0.4127716 sample estimates: mean of x mean of y 5.165973 6.033374

Chi-squared

Take a simple contingency table:

handed ------------ right left Total female RF LF male RM LM Total

For any values (RF,LF,RM,LM) we calculate the marginal totals and then we can ask, how surprising is the observed distribution given the marginal distributions?

Suppose the distribution is that we have 60 females and 80 right-handers, out of a total of 100 subjects. (Actually, the incidence of left-handers is about half) this much.

If hand bias and sex are not associated, what distribution of distributions do we expect by chance? Let's do a simulation in Python. Here are the results from 50,000 reps.

mean sd 95% CI RF 48.0 1.97 44.14 51.85 LF 12.0 1.97 8.15 15.86 RM 32.0 1.97 28.15 35.86 LM 8.0 1.97 4.14 11.85

I was quite surprised by this at first. On reflection, I guess it is implicit in the whole idea of the chi-square test, but it just completely passed me by that the standard deviation for each subgroup is necessarily the same.

It is connected to the fact that, as they say, there is only one degree of freedom. If our group of 60 females happens to have not 48 but 53 right-handers, then the rest is all determined. There will be 7 left-handed females, and for the 40 males there will be not 32 right-handers but only 27, and not 8 left-handers but 13. Each group will be shifted from its mean value by the same amount. So the total amount of "shifting" or variation from the mean, is the same for each group. And of course the sum of the square of this difference for the whole series of trials will also be the same. Huh!

import random,math random.seed(157) RF = list() LF = list() RM = list() LM = list() LL = [RF,LF,RM,LM] nameL = ['RF','LF','RM','LM'] def initOne(): R = 80; F = 60; N = 100 bias = list('R'*R + 'L'*(N-R)) sex = list('F'*F + 'M'*(N-F)) random.shuffle(sex) random.shuffle(bias) return bias,sex for i in range(50000): bias,sex = initOne() L = list() while sex: L.append(bias.pop() + sex.pop()) for i,case in enumerate(nameL): LL[i].append(L.count(case)) for i,L in enumerate(LL): #print L print nameL[i] + ' ', n = len(L) m = sum(L)*1.0/n sumsq = sum([(e-m)**2 for e in L]) sd = math.sqrt(sumsq/n) print str(m).rjust(4), print str(round(sd,2)).rjust(7), print str(round(m-1.96*sd,2)).rjust(7), print str(round(m+1.96*sd,2)).rjust(7)

Approaching Normal

This post continues with the theme of exploring the normal distribution. I wanted to take a look at how fast repeated samples from a uniform distribution approach the normal, and in the process exercise my R skills a bit. If you have not used R, I recommend it for graphics and exploring statistics. It works well (and looks beautiful) on my Mac. In the course I taught in Spring, there was a student who had trouble installing it under Windows, but that machine was hosed anyway. Disclaimer: I am not an R guru. If I'm not "doing it right", let me know.

The first thing is to model rolling a standard six-sided die. We want to sample from the uniform distribution of integers between 1 and 6. Don't forget to sample with replacement. As usual, I use yellow background for my code, and blue background to show what the program prints to the screen.

d = 1:6 u = sample(d,10000,replace=T) mean(u) var(u)

As we expect, the mean or expected value is 1/6 * (1 + 2 + 3 + 4 + 5 + 6) = 21/6 = 3.5 and the variance is 1/6 * 2 * (0.5² + 1.5² + 2.5²) = 1/3 * (0.25 + 2.25 + 6.25) = 1/3 * 8.75 = 2.92, and the (population) standard deviation is √2.92 = 1.71.

> mean(u) [1] 3.4786 > var(u) [1] 2.917434 > sd(u) [1] 1.708050 > summary(u) Min. 1st Qu. Median Mean 3rd Qu. Max. 1.000 2.000 3.000 3.479 5.000 6.000

We will take a look at the distribution of the numbers using the hist function. There are a couple of details to consider when using hist. (I often have to remind myself about its arguments by doing "?hist"). One is the argument "breaks"

breaks one of: a vector giving the breakpoints between histogram cells, a single number giving the number of cells for the histogram, a character string naming an algorithm to compute the number of cells (see ‘Details’), a function to compute the number of cells. In the last three cases the number is a suggestion only.

I often specify the number of cells explicitly, especially when I want to compare multiple plots. Since I had a little trouble with this plot, I tried a couple of different things, and I want to plot them all in the same window. For this, I use the formatting command "par" and tell it to make a set of plots in 1 row and 3 columns.

par(mfrow=c(1,3)) hist(u,breaks=6,col='blue') hist(u,right=F, breaks=6,col='blue') hist(u-0.5,breaks=6,col='blue')

As you can see, the first two histograms look funny. What is going on is that R is trying to bin the numbers in the vector with breakpoints exactly on the integer values 1, 2, 3... Since the numbers are themselves drawn from 1, 2, 3... it has to go right or left, and at the boundaries of the plot it looks weird. One argument controlling how this works is:

right logical; if TRUE, the histograms cells are right-closed (left open) intervals.

The default setting is "TRUE", which has cells formed as "right-closed (left-open) intervals"---whatever that means. The result is that at the left boundary all the values "1" and "2" have been binned together. Neither setting for "right" is what we want. My solution was to shift all the values to the left by 0.5 and then plot the result. The "1"s are plotted in the cell between 0 and 1.

Now, let's see what happens if we roll the dice again. What we're going to do is start with a vector of the size we need called z, that is filled with zeroes. (It has a length of 50,000). We need to initialize it because we will be updating at each round. We roll the dice six times and plot the results at each stage.



As you can see, we already have a reasonable approximation to the normal after summing just two numbers drawn at random from the "sample space" of 1, 2... , 6. And by n=6 the approximation is very good. Since the standard deviation goes as √variance, and the variances add, by n=6 we have a range of 31 (from 6 to 36) but the sd is only √(6*2.91) = 4.18 (4.19 in the figure).

Here is the code:

d = 1:6 par(mfrow=c(3,2)) z = rep(0,50000) for (i in 1:6) { B = i*6 z = z + sample(d,50000,replace=T) w = z+0.5 hist(w,breaks=B, xlim=c(min(w)-2,max(w)+2), col='gray90',freq=F, main=paste('sd = ', round(sd(w),2),sep=""), xlab=paste('round',i)) plot(function(x) dnorm(x,mean(w),sd(w)), 0,max(w),lwd=2,add=T,col='red') }

I found the code for the function that plots the normal distribution with the same mean and sd as our samples somewhere on the web. It is doing something a bit funny. I think it is making an "anonymous" (i.e. unnamed) function to feed to plot, and then applying that function with bounds = 0 and max(w), but I'm not real clear about how this works. The parameter lwd is for line width.

Finally, let's look a little more closely at the z vector after 6 rounds. We plot it in the same histogram as the normal density of the same mean and sd, switching which one is in back as we move from the left panel to the right.



The code:

mean(z) sd(z) summary(z) x=rnorm(length(z),mean(z),sd(z)) summary(x) par(mfrow=c(1,2)) hist(z,breaks=40,freq=F, xlim=c(0,40)) hist(x,col='gray80',breaks=40, freq=F,add=T) hist(x,col='gray80',breaks=40, freq=F,xlim=c(0,40)) hist(z,col='white', breaks=40,freq=F,add=T)

Notice, however, that the correspondence is not perfect, as shown by the results of calling the summary function.

> mean(z) [1] 21.00562 > sd(z) [1] 4.155018 > summary(z) Min. 1st Qu. Median Mean 3rd Qu. Max. 7.00 18.00 21.00 21.01 24.00 36.00 > x=rnorm(length(z),mean(z),sd(z)) > summary(x) Min. 1st Qu. Median Mean 3rd Qu. Max. 3.297 18.210 20.980 20.990 23.780 38.350

Our vector z seems a little fatter at the first and third quartiles and yet its minimum and maximum values are closer to the mean than the true normal distribution.

Normal approximation to the binomial

I know that the normal can be used as an approximation to the binomial. I was looking for a derivation of this, and I found it via google in a math forum. Doctor Anthony begins:

Derivation of the Normal distribution from the Binomial distribution --------------------------------------------------------------------- Let a variate take values 0, k, 2k, 3k, ..., nk with probabilities given by successive terms of (q + p)^n.

What's with the k? Well, we're eventually going to want non-integer terms. The expansion of (q + p)ⁿ is familiar:

qn + nqn-1 p + n(n-1)/2 qn-2 p2 + ...

The i^th term of the expansion is C(n,i).

Then the mean m = npk and the variance s^2 = npqk^2

OK. Notice use of the multiplication rule for variance from the other day.

Suppose: y = probability of occurrence of rk = C(n,r) p^r q^(n-r) Also let: y' = probability of occurrence (r+1)k = C(n,r+1)p^(r+1) q^(n-r-1) Then: y' - y = C(n,r+1)p^(r+1) q^(n-r-1) - C(n,r)p^r q^(n-r) n!p^r q^(n-r-1) = ---------------[(n-r)p - (r+1)q] (r+1)! (n-r)!

Hmm... I know that

y = [n! / r! (n-r)!] pr qn-r y' = [n! / (r+1)! (n-r-1)!] pr+1 qn-r-1 y' - y = ?

Lucie, you got some factoring to do. Let's deal with q and p first.

The left term has q^n-r-1 and the right term has q^n-r, so we can factor out
q^n-r-1, leaving a factor of q on the right-hand term in the brackets.

Similarly we can factor out p^r from both sides leaving a factor of p on the left.

The combination expressions expand as shown above. We can factor out n! from both sides. We can factor out 1/(r+1)! from both sides, if we first multiply top and bottom of the right-hand term by (r+1), leaving (r+1) on the top.

Similarly, we can factor out (n-r)! from both sides, if we first multiply top and bottom of the left-hand term by (n-r), leaving (n-r) behind on the top. So everything checks out so far. Next, he wants to divide by y:

And: y' - y 1 1 ------ = ------[np - r(p+q) - q] = ------[np - r - q] y (r+1)q (r+1)q (Equation 1)

Hmm...again. We're dividing the expression we had above by y.

y = [n! / r! (n-r)!] pr qn-r

We have:

n!pr qn-r-1 y' - y = ---------------[(n-r)p - (r+1)q] (r+1)! (n-r)!

So both n! and (n-r)! terms cancel. We also cancel r!, leaving a factor of (r+1) on the bottom. The p^r cancels, and the q^n-r also cancels leaves a factor of q on the bottom. So I get:

y' - y 1 ------ = ------[(n-r)p - (r+1)q] y (r+1)q

Now we have to figure out how to rearrange the term in brackets:

[(n-r)p - (r+1)q]

Expand, and then substitute for p + q = 1:

np - rp - rq - q np - r(p+q) - q np - r - q

It checks out. Doctor Anthony continues:

Let: x = rk - npk, so that x is now the variate measured from the mean. Then: r = x/k + np and r+1 = x/k + np + 1 Thus: k(r+1) = x + k + npk k^2 (r+1)q = (x + k + npk)qk

So far so good.

Multiply top and bottom of the righthand side of Equation 1 by k^2. Then: y' - y [(np-r)k - qk]k ------ = --------------- [note that (np-r)k = -x] y [x + k + npk]qk

Go back to what we had, and then multiply top and bottom by k²:

y' - y [np - r - q] k^2 ------ = ---------------- y (r+1)q k^2

Hmm... The top is fine, but on the bottom we had

(r+1) q k2

We need to get to:

[x + k + npk]qk

He says:

[note that (np-r)k = -x]

OK, so we have:

(r+1) q k2 (rk + k) q k Since: (np - r) k = -x rk = npk + x Substituting: (x + k + npk) q k

Moving on to substitute for (np-r) k = -x on top and multiplying out on the bottom yields:

(-x - kq)k = ---------------- npqk2 + (x+k)qk

Finally, we now let k = dx, so that y' - y = dy and let n ->infinity in such a way that nk^2 is finite. Equation 2 can then be written as: dy (-x - q dx)dx ---- = ---------------- y s^2 + (x+dx)q dx

The only tricky part here was that we've replaced npqk² by s².
Now he says:

As dx -> 0 this becomes: dy -x dx ---- = ------ y s^2

And we're there! If we integrate the left side we get ln(y), and the right side is
-x² / 2s²

y = A exp { -x² / 2s² }

Gaussian normalizing constant

Last time we followed a derivation for the normal or Gaussian distribution that gave us the general form:

p(x) = A exp { -x2/2V }

What we need to do now is to evaluate the constant A. As I'm sure you remember, it has something to do with π, more accurately with √(2π). The constant doesn't really change the distribution (e.g. the fraction of the density that lies between x = +σ and x = +∞). What the constant does is to "normalize" the values so that the total density over the range of x (from -∞ to +∞) to equals one, a requirement for a real probability density function. It squishes the graph in the y-dimension. It is probably obvious that, if the total value (∫ exp { -x²/2V } over the whole range = 1/A, then multiplying by A will normalize it.

The V in the formula above turns out to be the variance of the distribution. It's usually written as σ² but I want to leave it as it is. Let us start with V = 1 and then explore the influence of V separately. That means σ also equals one.

We can evaluate the constant A easily using R. Since the function is symmetrical, we can integrate between 0 and +∞, and then multiply the result by two. We're going to do two things to simplify our lives. First, rather than integrate, we will sum a bunch of small intervals. Second, rather than go to +∞, we will only go as far as there is significant density. Remembering the rules for z-scores, you can guess that we'll get 99% accuracy if we go out to 3*σ. Let's aim even higher and go to 5*σ.

T = 5 I = 0.1 x = seq(0,T,by=I) x = x + I/2 x = x[-length(x)] head(x)

> head(x) [1] 0.05 0.15 0.25 0.35 0.45 0.55

We specify the intervals using the variable I. Here I = 0.1, but you can make it smaller if you wish. We construct a vector v which has the values between 0 and T with a spacing of I = 0.1. The sums are more accurate if we bump out all the values in the vector by I/2 (and then remove the last one, which is outside the interval we want to measure. A single interval looks something like this.




We approximate the area of the triangle on the left (red side) by the area of the triangle on the right. We calculate the vector of areas of all the slices and store the result in y, then sum over y. We multiply by 2 to adjust for the symmetrical values for -x.

y = I * exp(-0.5*x**2) k = 2*sum(y) k sqrt(2*pi)

Because we know the true value, we can compare our answer it.

> k [1] 2.506627 > sqrt(2*pi) [1] 2.506628

We can also generalize the above code by explicitly including V. Trying different values of V, it is clear that the constant k is actually √(2 π V), and the A in the density function is 1/k. Note the use of √V = σ in the definition of T.

V = 16 T = 5*sqrt(V) I = 0.1 x = seq(0,T,by=I) x = x + I/2 x = x[-length(x)] y = I * exp(-0.5*x**2/V) k = 2*sum(y) sqrt(2*pi) sqrt(2*pi*V)

> k [1] 10.02651 > sqrt(2*pi) [1] 2.506628 > sqrt(2*pi*V) [1] 10.02651

Dan Teague's derivation (pdf) evaluates A using calculus. I just wanted to indicate the way π creeps in. Here are some equations from him:



The double integral over x and y (equation 3) is converted to an integral in radial coordinates. The ∫ θ then contributes the π And since the conversion to radial coordinates involved x times y, going back to p(x) gives us √π. I don't know enough to be sure this is OK, but it certainly looks reasonable.

The normal (Gaussian) distribution

I'm not very good at proofs, but I wanted to try to understand where the normal distribution comes from. In fact, we saw in an earlier post that we can show by simulation that the Central Limit Theorem seems to be correct. Regardless of the underlying distribution, the sample mean x is normally distributed if the sample size is sufficiently large.

However, let's try this argument, which is originally due to Sir John F. W. Herschel.

Imagine that you are throwing darts at the origin of the x,y plane. Under perfect conditions, you would hit the center dead on every time. However, conditions aren't perfect. The wind is gusting, the music is loud, your blood alchohol is modestly elevated, there are other distractions. As a result, small errors creep in and the pattern over time looks like so:



The R code:

x=rnorm(1000) y=rnorm(1000) L=c(-3,3) plot(x,y,pch=16,xlim=L,ylim=L,col='blue') lines(c(0,0),c(-3,3),lty=2,lwd=2) lines(c(-3,3),c(0,0),lty=2,lwd=2)

Now, there is some unknown function for the probability that a dart will land in the interval between x and x + ∆x. Obviously, the probability depends on x, with a maximum at x = 0 and then decreasing to zero as x gets large. We designate that function as a probability density function p(x) and evaluate the density over the interval to get the probability that the dart lands in the interval:

Prob = p(x) ∆x

Now we consider a small area of size ∆x∆y. If:

the errors in perpendicular directions are independent

then we expect that p(x) = p(y) and we can get the probability that a dart lands in the small rectangle bounded by x, y and x + ∆x, y + ∆y as:

Prob = p(x)∆x p(y)∆y

In fact, if we assume that the errors do not depend on the orientation of the coordinate system, then the probability is a function only of r, the radial distance from the origin, so we can write

Prob = g(r)∆x∆y g(r)∆x∆y = p(x)∆x p(y)∆y g(r) = p(x) p(y)

This assumption of rotational independence will lead us directly to the answer, as you will see. As Hamming says, since r does not depend on the angle θ, (but x and y do), we can take the partial derivative with respect to θ of g(r) and set it equal to zero, so that:



We can parse this. We used the standard multiplication rule (twice): "this times the derivative of that plus that times the derivative of this." We use it to generate the first line (taking the partial derivative of p(x) p(y)). And then, we need to actually find the partial derivatives of p(x) and p(y) with respect to θ, where x = r cos(θ) and y = r sin(θ). We use the multiplication rule again, and the fact that the derivative of the sine is just the cosine, while the derivative of the cosine is minus the sine. Thus, for example, the partial derivative of x with respect to θ is simply -y.

As stated, this gives:

p(x) p'(y)(x) - p(y)p'(x)(y) = 0 p'(x)/x p(x) = p'(y)/y p(y)

Since x and y are both variables:

p'(x)/x p(x) = p'(y)/y p(y) = K p'(x)/p(x) = Kx

We need a function p(x) whose derivative p'(x) is equal to p(x) times x times a constant. Remember the exponential function from a few days ago?

p(x) = A exp { Kx2/2 }

Since we assume that large errors are less likely than small ones, K < 0, so we can define another constant V = - 1/K and

p(x) = A exp { -x2/2V }

This is the normal distribution with variance V.

It is amazing how far we got with this argument! We assumed:

(1) the errors do not depend on the orientation of the coordinate system.
(2) errors in perpendicular directions are independent. This means that being too high doesn't alter the probability of being off to the right.
(3) large errors are less likely than small errors.

The pdf from Dan Teague has more. Notice that although we started talking about a probability distribution in two dimensions, the function we end up with is for one dimension.

Even better, James Clerk Maxwell used the same argument in three dimensions to derive his expression for the distribution of molecular velocities in a gas. Here is a very cool simulation that shows the distribution.

Statistical doodling: variance

In Bolstad, Chapter 5, there is a proof of the following statement about the variance of independent random variables X and Y.

Var(X + Y) = Var(X) + Var(Y)

There is a lot more discussion here. The post calls this the "Pythagorean Theorem of Statistics", since an equivalent formulation is:

SD2(X + Y) = SD2(X) + SD2(Y)

I don't want to detail the proof, but I did fool around a bit in R to explore this:

set.seed(1357) u = rnorm(10000,5,1) var(u) var(u + 7) var(u-250) var(3*u) var(u/5)

Here is what it prints:

> var(u) [1] 0.9962947 > var(u + 7) [1] 0.9962947 > var(u-250) [1] 0.9962947 > var(3*u) [1] 8.966652 > var(u/5) [1] 0.03985179

So, if we add or subtract a constant C, the variance is unchanged. But if we multiply by C, the variance is multiplied by C²; and if we divide by C, the variance is divided by C².

Now consider a second set of numbers from rnorm. The first vector has a mean of 5 and sd of 2 (variance of 4), while the second has a mean of 4 and sd of 3 (variance of 9).

u = rnorm(1000,5,2) v = rnorm(1000,4,3) var(u+v) var(u-v)

> u = rnorm(1000,5,2) > v = rnorm(1000,4,3) > var(u) [1] 3.99337 > var(v) [1] 9.470766 > var(u+v) [1] 12.76349 > var(u-v) [1] 13.65514

Our simulation confirms the rule that the variances add.

And finally, look at multiplication:

u = rnorm(1000,0,1) v = rnorm(1000,0,1) var(u*v) var((u+1)*v) var((u+2)*v) var((u+3)*v) var((u+2)*(v+2))

The variance depends on the mean of the distributions. Here, the variances of u and v (as well as u + 1,2..3) are always 1. For means of:

0,0: var = 1.0 1,0: var = 2.2 2,0: var = 5.5 3,0: var = 10.9 2,2: var = 9.4

I found an expression here:

Var(XY) = Var(X)*Var(Y) + Var(X)*E[Y]^2 + E[X]^2 Var(Y)

That is:

v(X*Y) = vX*vY + vX*mY2 + mX2*vY

In the cases above (variance is unchanged and equal to 1) we have:

0,0: v(X*Y) = 1 + 1*0 + 0 *1 = 1 1,0: v(X*Y) = 1 + 1*1 + 0 *1 = 2 2,0: v(X*Y) = 1 + 1*22 + 0 *1 = 5 3,0: v(X*Y) = 1 + 1*32 + 0 *1 = 10 2,2: v(X*Y) = 1 + 1*22 + 22*1 = 9

Looks correct.

Sweet spots of the bell curve

Continuing with the theme of "basic stuff I never learned," here is something interesting about the normal distribution. It turns out that the inflection points of the bell curve are the points where x = σ. I think that's pretty amazing. Let's see if we can prove it using a small bit o'calculus.

As we come over the top of the curve and head down, the slope is becoming increasingly negative. But at some point the slope reaches its maximum negative value and then starts to turn less negative (more positive). At one instance the slope of the slope or second derivative of the pdf is zero. So we need to differentiate the normal density function twice, set it equal to zero, and then solve.

I have to admit I got too confused in the middle of the calculation, so I needed help. I googled 'second derivative normal distribution' and found this.



The pdf for the normal distribution is equation (1). We're going to differentiate twice and set that result equal to zero, so the constant out front can be ignored. We rewrite the pdf as equation (2). To simplify the notation, we will define f(x) as in equation (3) and then we can rewrite equation (2) as equation (4).

We will use the result in equation (5) several times. This is just a generalization of what I mentioned the other day with respect to the exponential distribution and its cdf.

We consider the exponent part as f(x) in (3) and (4). We use the chain rule to find its derivative. Set x - μ / σ = g(x) and then do: df/dx = df/dg * dg/dx and we obtain equation (6).

We use the results from (5) and (6) in figuring out the derivative of φ(x) as shown in equation (7). As mentioned, the derivative of exp { f(x) } is just f '(x) exp { f(x) }. We calculated f'(x) in (6) and we put the results together in (7). Now we have something substantially more complicated but it is really just the product of two functions each of which we know how to differentiate.

Remember that the derivative of g(x) f(x) is g'(x) f(x) + g(x) f '(x). "This times the derivative of that plus that times the derivative of this."



In the second panel we restate (7). To do the second differentiation, we first pull the constant out front and do this (x-μ / σ) times the derivative of that (7), plus that (exp { f(x) }) times the derivative of this (1/sigma;). We pull out the common factors and obtain (9).

We set this equal to zero, and now it is easy to see that the solutions are as shown in (11). Pretty neat!

And if anybody can tell me how to make my beautiful images not look like crap in blogger I would appreciate it. I suppose I need to RTFM.

Exponential density 3

In example 4.20 of Grinstead and Snell there is a nice conjunction of Bayes theorem and use of the exponential density. Recall that the exponential pdf is:



To find the probability that X happens (hard-drive failure, radioactive decay) within a certain time period, we integrate the pdf over the interval. For example, the probability that the failure happens after a particular time t is:



(Naturally, since the cdf(t) is 1 minus this value).

Here is one form of Bayes theorem:



Now, consider two events E and F defined as follows:

E is the event that failure happens after time r F is the event that failure happens after time r + s

Note that P(F and E) = P(F) because F is totally contained within E.

Then:



The probability of failure after time r + s, when we know that failure occurs after time r, does not depend on r at all but is only a function of s. This is the memoryless property of the exponential function, alluded to in a previous post.

Exponential density 2

Think of e as a function (rather than the irrational number 2.71828...).

In R:

e <- function(x) { 2.71828**(x) } plot(e,0,5,lwd=10,col='gray70') plot(exp,0,5,lwd=2,col='red',add=T)

Then, e can be defined as the function whose derivative is itself. To see this:

plot(exp,0,5,lwd=5,xlim=c(0,5),ylim=c(0,75)) par(lwd=2) colors=c('blue','darkgreen','red','magenta') f <- function(x) { y=exp(x) points(x,y,col=colors[x],pch=16,cex=2) lines(c(x-1,x+1),c(y-y,y+y), col=colors[x]) } for (i in 1:4) f(i)

What did we do? First, we plotted the exponential function exp between 0 and 5 (heavy black line). Then, for each x in the series 1:4, we calculated y=exp(x) and plotted the point x,y in color. At the same time, we constructed a line with slope (y+y)/2 = y, passing through x,y. It is clear that the slope of the curve is equal to the value of the function at that point.

Another fun way to see this is to look at the infinite series for e^x:



Can you see that the derivative with respect to x of this series is identical to the series itself?

Using simple calculus (the chain rule), we can also show that



Looking at the probability density function (pdf) for the exponential distribution,



and cumulative distribution function, we confirm that the pdf is the derivative of the cdf, as it should be. The cdf is:

Exponential density

I've been reading (and re-reading) An Introduction to Probability by Grinstead and Snell. It is a wonderful book, available from here as a pdf. It goes slowly, has lots of explanation and many problems, as well as interesting historical perspective. I like it so much that I bought a hard copy.

I'm trying to understand the exponential density better. In example 2.17 of the book they pose the problem of modeling the time-to-breakdown of a hard drive by the exponential density:



If the average time-to-breakdown is 30 months, and we have already run the computer for 15 months with no breakdown, what is the current expected time-to-breakdown?

We use R to explore the question. The R function rexp gives random samples from the exponential density with a rate parameter r (the inverse of lambda above). Think of each of these as a possible lifetime for our drive. Since we know that the lifetime exceeds 15 months, filter the vector x for values > 15 and save in y.

r = 1/30 x = rexp(100000,rate=r) sel = x > 15 y = x[sel]

Plot histograms of the density (not counts, which R calls freq).

hist(y,breaks=100,xlim=c(0,150),freq=F) hist(x,breaks=100,col='gray70',freq=F,add=T)

It is clear that the y distribution is the same as x, just shifted over by 15.



We confirm this by looking at summary statistics (adjusting y first by subtracting 15). We see they are essentially identical:

> summary(x) Min. 1st Qu. Median 4.915e-05 8.580e+00 2.075e+01 Mean 3rd Qu. Max. 3.001e+01 4.170e+01 3.780e+02 > summary(y-15) Min. 1st Qu. Median 1.309e-05 8.627e+00 2.087e+01 Mean 3rd Qu. Max. 3.004e+01 4.185e+01 3.630e+02 freq=F,add=T)