Login
📚 Statistical Inference for Everyone
Chapters ▾
1 Introduction to Probability2 Applications of Probability3 Random Sequences and Visualization4 Introduction to Model Comparison5 Applications of Model Comparison6 Introduction to Parameter Estimation7 Priors, Likelihoods, and Posteriors8 Common Statistical Significance Tests9 Applications of Parameter Estimation and Inference10 Multi-parameter Models11 Introduction to MCMC12 Concluding Thoughts13 Appendix A: Computational Analysis14 Appendix B: Notation and Standards15 Appendix C: Common Distributions and Their Properties16 Appendix D: Tables

9.5 Computer Examples

from sie import *

Iris Example

data=load_data('data/iris.csv')
x_sertosa=data[data['class']=='Iris-setosa']['petal length [cm]']
x_virginica=data[data['class']=='Iris-virginica']['petal length [cm]']
x_versicolor=data[data['class']=='Iris-versicolor']['petal length [cm]']
print x_sertosa[:10]  # print the first 10
Show expected output
0    1.4
1    1.4
2    1.3
3    1.5
4    1.4
5    1.7
6    1.4
7    1.5
8    1.4
9    1.5
Name: petal length [cm], dtype: float64
x=x_sertosa
mu=sample_mean(x)
N=len(x)
sigma=sample_deviation(x)/sqrt(N)
t_sertosa=tdist(N,mu,sigma)

print "total number of data points:",N
print "best estimate:",mu
print "uncertainty:",sigma
Show expected output
total number of data points: 50
best estimate: 1.464
uncertainty: 0.0245381834898
x=x_versicolor
mu=sample_mean(x)
N=len(x)
sigma=sample_deviation(x)/sqrt(N)
t_versicolor=tdist(N,mu,sigma)

print "total number of data points:",N
print "best estimate:",mu
print "uncertainty:",sigma
Show expected output
total number of data points: 50
best estimate: 4.26
uncertainty: 0.0664554477121
x=x_virginica
mu=sample_mean(x)
N=len(x)
sigma=sample_deviation(x)/sqrt(N)
t_virginica=tdist(N,mu,sigma)

print "total number of data points:",N
print "best estimate:",mu
print "uncertainty:",sigma
Show expected output
total number of data points: 50
best estimate: 5.552
uncertainty: 0.078049696361
distplot2([t_sertosa,t_versicolor,t_virginica],show_quartiles=False)
Show expected output
<matplotlib.figure.Figure at 0x1058d9690>
Three Student-t posterior curves for iris petal lengths from the sample code: a tall narrow blue peak near 1.45 cm (Setosa), and two broader peaks near 4.25 cm (green) and 5.55 cm (red).
distplot(t_virginica)
Posterior for the Virginica mean petal length: a bell curve centered at 5.55 cm with dashed percentile lines from 5.36 (1%) to 5.74 (99%).
credible_interval(t_versicolor)
Show expected output
(4.1265203051077082, 4.2599999999999998, 4.3934796948922914)
credible_interval(t_virginica)
Show expected output
(5.3952325713636533, 5.5519999999999996, 5.7087674286363459)

Sunrise

dist=beta(h=365,N=365)
distplot(dist)
Steeply rising posterior for the probability that Virginica petal length exceeds Versicolor's, crowded against 1.0: the 1% percentile is already 0.99 and all higher percentiles round to 1.00.
credible_interval(dist)
Show expected output
(0.98997171634278669, 0.99810794743679487, 0.99993082805373457)

Cancer Example

dist=beta(h=7,N=10)
distplot(dist,figsize=(8,5))
Posterior for the proportion θ from the sample code: a left-skewed hump peaking near 0.7 with dashed percentile lines running from 0.34 (1%) to 0.92 (99%).
credible_interval(dist)
Show expected output
(0.39025744042757882, 0.67619553741481253, 0.89073655618090186)

Essentially no evidence of any effect over 50 percent.

Pennies

data1=load_data('data/pennies1.csv')
print data1
year,mass=data1['Year'],data1['Mass [g]']
Show expected output
    Year  Mass [g]
0   1960     3.133
1   1961     3.083
2   1962     3.175
3   1963     3.120
4   1964     3.100
5   1965     3.060
6   1966     3.100
7   1967     3.100
8   1968     3.073
9   1969     3.076
10  1970     3.100
11  1971     3.110
12  1972     3.080
13  1973     3.100
14  1974     3.093
plot(year,mass,'o')
xlabel('year')
ylabel('Mass per Penny [g]')
Show expected output
<matplotlib.text.Text at 0x1087c2d90>
Scatter of the mass of US pennies by year, 1960-1974: points fluctuate between about 3.06 and 3.18 grams with no visible trend; most years sit near 3.10 grams.
x=mass
mu=sample_mean(x)
N=len(x)
sigma=sample_deviation(x)/sqrt(N)
t_penny1=tdist(N,mu,sigma)

distplot(t_penny1,label='mass [g]')
Posterior distribution for the true penny mass: a bell curve peaking at 3.10 grams with dashed percentile lines from 3.08 (1%) to 3.12 (99%).
CI=credible_interval(t_penny1,percentage=99)
print CI
Show expected output
(3.0790129206702002, 3.1002000000000001, 3.1213870793298)
plot(year,mass,'o')
credible_interval_plot(t_penny1,percentage=99)
xlabel('year')
ylabel('Mass per Penny [g]')
Show expected output
<matplotlib.text.Text at 0x1087fcf10>
The 1960-1974 penny-mass scatter with the best estimate drawn as a dashed horizontal line at 3.10 grams and dotted uncertainty bands just above and below it; all but two points fall inside the bands.

Do the 2 datasets

data2=load_data('data/pennies2.csv')
print data2
year1,mass1=year,mass
year2,mass2=data2['Year'],data2['Mass [g]']
Show expected output
    Year  Mass [g]
0   1989     2.516
1   1990     2.500
2   1991     2.500
3   1992     2.500
4   1993     2.503
5   1994     2.500
6   1995     2.497
7   1996     2.500
8   1997     2.494
9   1998     2.512
10  1999     2.521
11  2000     2.499
12  2001     2.523
13  2002     2.518
14  2003     2.520
x=mass1
mu=sample_mean(x)
N=len(x)
sigma=sample_deviation(x)/sqrt(N)
t_penny1=tdist(N,mu,sigma)

x=mass2
mu=sample_mean(x)
N=len(x)
sigma=sample_deviation(x)/sqrt(N)
t_penny2=tdist(N,mu,sigma)

distplot2([t_penny1,t_penny2],show_quartiles=False,label='mass [g]')
legend([r'$\mu_1$',r'$\mu_2$'])
Show expected output
<matplotlib.figure.Figure at 0x1087d3f10>
Show expected output
<matplotlib.legend.Legend at 0x1088198d0>
Two sharply separated posteriors for the mean penny mass, labeled μ1 and μ2: the μ2 curve is a tall green spike at about 2.50 grams (post-1975 pennies) and the μ1 curve a shorter blue peak at about 3.10 grams (earlier pennies).
plot(year1,mass1,'o')
credible_interval_plot(t_penny1,percentage=99)
plot(year2,mass2,'ro')
credible_interval_plot(t_penny2,percentage=99,xlim=[1989,2005])
xlabel('year')
ylabel('Mass per Penny [g]')
Show expected output
<matplotlib.text.Text at 0x10907e310>
Penny masses from 1960 to 2005 in two clusters: blue points near 3.10 grams before 1975 and red points near 2.50 grams after 1989, each cluster with its own dashed best-estimate line and dotted uncertainty band. The two bands are far apart.

Distribution of the difference, normal approximation

N1=len(mass1)
N2=len(mass2)

mu1=sample_mean(mass1)
mu2=sample_mean(mass2)


sigma1=(1+20.0/N1**2)*sample_deviation(mass1)/sqrt(N1)
sigma2=(1+20.0/N2**2)*sample_deviation(mass2)/sqrt(N1)


delta_12=mu1-mu2
sigma_delta12=sqrt(sigma1**2+sigma2**2)

dist_delta=normal(delta_12,sigma_delta12)
distplot(dist_delta)
Posterior for a proportion from the sample code: a bell curve centered at 0.59 with dashed percentile lines from 0.57 (1%) to 0.61 (99%).

clearly larger than zero at well over the 99

Ball Bearing Sizes

data1=[1.18,1.42,0.69,0.88,1.62,1.09,1.53,1.02,1.19,1.32]
data2=[1.72,1.62,1.69,0.79,1.79,0.77,1.44,1.29,1.96,0.99]
N1=len(data1)
N2=len(data2)
mu1=sample_mean(data1)
mu2=sample_mean(data2)
print mu1,mu2
Show expected output
1.194 1.406
S1=sample_deviation(data1)
S2=sample_deviation(data2)
print S1,S2
Show expected output
0.289681817786 0.428309337849
sigma1=S1/sqrt(N1)
sigma2=S2/sqrt(N2)
print sigma1,sigma2
Show expected output
0.091605434094 0.135443305072
dist1=normal(mu1,sigma1)
dist2=normal(mu2,sigma2)
distplot2([dist1,dist2],show_quartiles=False,label='size [microns]')
legend([r'$\mu_1$',r'$\mu_2$'])
Show expected output
<matplotlib.figure.Figure at 0x105d61390>
Show expected output
<matplotlib.legend.Legend at 0x108ca1ad0>
Two overlapping posteriors for ball-bearing sizes labeled μ1 and μ2: a taller blue curve peaking at about 1.2 microns and a broader green curve peaking at about 1.4 microns, with substantial overlap between them.

Adapted from Statistical Inference for Everyone, by Brian Blais (Bryant University), licensed under CC BY-SA 4.0 (dual-licensed under the GNU FDL 1.2 or later; this adaptation uses the CC BY-SA grant). Changes were made; this adaptation is distributed under the same license. License: CC-BY-SA-4.0.