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.