Skip to content

Commit a77a7b8

Browse files
ericmjlericmjl
committed
first commit of bonus notebook on finches
1 parent 99c4f3d commit a77a7b8

2 files changed

Lines changed: 194 additions & 1 deletion

File tree

‎notebooks/bonus-exploration-finches.ipynb‎

Lines changed: 191 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,191 @@
1+
{
2+
"cells": [
3+
{
4+
"cell_type": "code",
5+
"execution_count": null,
6+
"metadata": {},
7+
"outputs": [],
8+
"source": [
9+
"import pymc3 as pm\n",
10+
"import matplotlib.pyplot as plt\n",
11+
"import numpy as np\n",
12+
"from data import load_finches_2012, load_finches_1975\n",
13+
"from utils import ECDF\n",
14+
"\n",
15+
"%load_ext autoreload\n",
16+
"%autoreload 2\n",
17+
"%matplotlib inline\n",
18+
"%config InlineBackend.figure_format = 'retina'"
19+
]
20+
},
21+
{
22+
"cell_type": "code",
23+
"execution_count": null,
24+
"metadata": {},
25+
"outputs": [],
26+
"source": [
27+
"df12 = load_finches_2012()\n",
28+
"df12['shape'] = df12['beak_depth'] / df12['beak_length']\n",
29+
"\n",
30+
"df12 = df12[df12['species'] != 'unknown']\n",
31+
"df75 = load_finches_1975()\n",
32+
"\n",
33+
"df = df12 # convenient alias"
34+
]
35+
},
36+
{
37+
"cell_type": "code",
38+
"execution_count": null,
39+
"metadata": {},
40+
"outputs": [],
41+
"source": [
42+
"df12.head(5)"
43+
]
44+
},
45+
{
46+
"cell_type": "code",
47+
"execution_count": null,
48+
"metadata": {},
49+
"outputs": [],
50+
"source": [
51+
"fortis_idx = df[df['species'] == 'fortis'].index\n",
52+
"scandens_idx = df[df['species'] == 'scandens'].index"
53+
]
54+
},
55+
{
56+
"cell_type": "code",
57+
"execution_count": null,
58+
"metadata": {},
59+
"outputs": [],
60+
"source": [
61+
"# Mega-model incorporating shape as well. \n",
62+
"# We will also analyze the SD in addition to the mean.\n",
63+
"\n",
64+
"with pm.Model() as beak_model:\n",
65+
" # SD can only be positive, therefore it is reasonable to constrain to >0\n",
66+
" # Likewise for betas.\n",
67+
" sd_hyper = pm.HalfCauchy('sd_hyper', beta=100, shape=(2,))\n",
68+
" beta_hyper = pm.HalfCauchy('beta_hyper', beta=100, shape=(2,))\n",
69+
" \n",
70+
" # Beaks cannot be of \"negative\" mean, therefore, HalfNormal is \n",
71+
" # a reasonable, constrained prior.\n",
72+
" mean_depth = pm.HalfNormal('mean_depth', sd=sd_hyper[0], shape=(2,))\n",
73+
" sd_depth = pm.HalfCauchy('sd_depth', beta=beta_hyper[0], shape=(2,))\n",
74+
" \n",
75+
" mean_length = pm.HalfNormal('mean_length', sd=sd_hyper[1], shape=(2,))\n",
76+
" sd_length = pm.HalfCauchy('sd_length', beta=beta_hyper[1], shape=(2,))\n",
77+
"\n",
78+
" nu = pm.Exponential('nu', lam=1/29.) + 1\n",
79+
" \n",
80+
" # Define the likelihood distribution for the data.\n",
81+
" depth = pm.StudentT('depth', \n",
82+
" nu=nu,\n",
83+
" mu=mean_depth[df['species_enc']], \n",
84+
" sd=sd_depth[df['species_enc']], \n",
85+
" observed=df['beak_depth'])\n",
86+
" \n",
87+
" length = pm.StudentT('length',\n",
88+
" nu=nu,\n",
89+
" mu=mean_length[df['species_enc']],\n",
90+
" sd=sd_length[df['species_enc']],\n",
91+
" observed=df['beak_length'])\n",
92+
" \n",
93+
" shape = pm.Deterministic('shape', depth / length)"
94+
]
95+
},
96+
{
97+
"cell_type": "code",
98+
"execution_count": null,
99+
"metadata": {},
100+
"outputs": [],
101+
"source": [
102+
"with beak_model:\n",
103+
" trace = pm.sample(2000)"
104+
]
105+
},
106+
{
107+
"cell_type": "code",
108+
"execution_count": null,
109+
"metadata": {},
110+
"outputs": [],
111+
"source": [
112+
"pm.traceplot(trace, varnames=['mean_length', 'mean_depth'])"
113+
]
114+
},
115+
{
116+
"cell_type": "code",
117+
"execution_count": null,
118+
"metadata": {},
119+
"outputs": [],
120+
"source": [
121+
"pm.traceplot(trace, varnames=['sd_length', 'sd_depth'])"
122+
]
123+
},
124+
{
125+
"cell_type": "code",
126+
"execution_count": null,
127+
"metadata": {},
128+
"outputs": [],
129+
"source": [
130+
"samples = pm.sample_ppc(trace, model=beak_model)\n",
131+
"samples"
132+
]
133+
},
134+
{
135+
"cell_type": "code",
136+
"execution_count": null,
137+
"metadata": {},
138+
"outputs": [],
139+
"source": [
140+
"fig = plt.figure()\n",
141+
"ax = fig.add_subplot(111)\n",
142+
"x, y = ECDF((samples['depth'][:, fortis_idx] / samples['length'][:, fortis_idx]).flatten())\n",
143+
"ax.plot(x, y)\n",
144+
"x, y = ECDF(df.loc[fortis_idx, 'shape'])\n",
145+
"ax.plot(x, y)"
146+
]
147+
},
148+
{
149+
"cell_type": "code",
150+
"execution_count": null,
151+
"metadata": {},
152+
"outputs": [],
153+
"source": [
154+
"fig = plt.figure()\n",
155+
"ax = fig.add_subplot(111)\n",
156+
"x, y = ECDF(df['shape'])\n",
157+
"ax.plot(x, y, label='data')\n",
158+
"# x, y = ECDF(trace['shape'][0, :])\n",
159+
"# ax.plot(x, y, label='posterior')\n"
160+
]
161+
},
162+
{
163+
"cell_type": "code",
164+
"execution_count": null,
165+
"metadata": {},
166+
"outputs": [],
167+
"source": []
168+
}
169+
],
170+
"metadata": {
171+
"kernelspec": {
172+
"display_name": "bayesian-modelling-tutorial",
173+
"language": "python",
174+
"name": "bayesian-modelling-tutorial"
175+
},
176+
"language_info": {
177+
"codemirror_mode": {
178+
"name": "ipython",
179+
"version": 3
180+
},
181+
"file_extension": ".py",
182+
"mimetype": "text/x-python",
183+
"name": "python",
184+
"nbconvert_exporter": "python",
185+
"pygments_lexer": "ipython3",
186+
"version": "3.6.6"
187+
}
188+
},
189+
"nbformat": 4,
190+
"nbformat_minor": 2
191+
}

‎notebooks/data.py‎

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -11,7 +11,9 @@ def load_finches_2012():
1111

1212
def load_finches_1975():
1313
path = '../data/finch_beaks_1975.csv'
14-
return load_finches(path)
14+
df = load_finches(path)
15+
df = df.rename_column('beak_length_mm', 'beak_length').rename_column('beak_depth_mm', 'beak_depth')
16+
return df
1517

1618

1719
def load_finches(path):

0 commit comments

Comments
 (0)