forked from aloctavodia/Doing_bayesian_data_analysis
-
Notifications
You must be signed in to change notification settings - Fork 0
/
plot_post.py
70 lines (62 loc) · 2.99 KB
/
plot_post.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
from __future__ import division
import numpy as np
from scipy import stats
import matplotlib.pyplot as plt
from hpd import hpd
def plot_post(param_sample_vec, cred_mass=0.95, comp_val=False,
ROPE=False, ylab='', xlab='parameter', fontsize=14, labelsize=14,
title='', framealpha=1, facecolor='skyblue', edgecolor='white',
show_mode=True, bins=50):
#compute HDI
HDI = hpd(param_sample_vec, 1-cred_mass)
post_summary = {'mean':0,'median':0,'mode':0, 'hdi_mass':0,'hdi_low':0,
'hdi_high':0, 'comp_val':0, 'pc_gt_comp_val':0, 'ROPE_low':0,
'ROPE_high':0, 'pc_in_ROPE':0}
post_summary['mean'] = np.mean(param_sample_vec)
post_summary['median'] = np.median(param_sample_vec)
post_summary['mode'] = stats.mode(param_sample_vec)[0]
post_summary['hdi_mass'] = cred_mass
post_summary['hdi_low'] = HDI[0]
post_summary['hdi_high'] = HDI[1]
# Plot histogram.
n, bins, patches = plt.hist(param_sample_vec, normed=True, bins=bins,
edgecolor=edgecolor, facecolor=facecolor)
plt.xlabel(xlab, fontsize=fontsize)
plt.ylabel(ylab, fontsize=fontsize)
plt.title(title, fontsize=fontsize)
cv_ht = 0.75*np.max(n)
cen_tend_ht = 0.9 * cv_ht
ROPE_text_ht = 0.55 * cv_ht
# # Display mean or mode:
if show_mode:
plt.plot(0, label='mode = %.2f' % post_summary['mode'], alpha=0)
else:
plt.plot(0, label='mean = %.2f' % post_summary['mean'], alpha=0)
# Display the comparison value.
if comp_val is not False:
pc_gt_comp_val = 100 * np.sum(param_sample_vec > comp_val)/len(param_sample_vec)
pc_lt_comp_val = 100 - pc_gt_comp_val
plt.plot([comp_val, comp_val], [0, cv_ht], color='darkgreen',
linestyle='--', linewidth=2,
label='%.1f%% <%.1f < %.1f%%'
% (pc_lt_comp_val, comp_val, pc_gt_comp_val))
post_summary['comp_val'] = comp_val
post_summary['pc_gt_comp_val'] = pc_gt_comp_val
# # Display the ROPE.
if ROPE is not False:
rope_col = 'darkred'
pc_in_ROPE = round(np.sum((param_sample_vec > ROPE[0]) & (param_sample_vec < ROPE[1]))/len(param_sample_vec)*100)
plt.plot([ROPE[0], ROPE[0]], [0, 0.96*ROPE_text_ht], color=rope_col,
linestyle=':', linewidth=4,
label='%.1f%% in ROPE' % pc_in_ROPE)
plt.plot([ROPE[1], ROPE[1]], [0, 0.96*ROPE_text_ht], color=rope_col,
linestyle=':', linewidth=4)
post_summary['ROPE_low'] = ROPE[0]
post_summary['ROPE_high'] = ROPE[1]
post_summary['pc_in_ROPE'] = pc_in_ROPE
# # Display the HDI.
plt.plot(HDI, [0, 0], linewidth=6, color='k', label='HDI %.1f%% %.3f-%.3f' % (cred_mass*100, HDI[0], HDI[1]))
plt.legend(loc='upper left', fontsize=labelsize, framealpha=framealpha)
frame = plt.gca()
frame.axes.get_yaxis().set_ticks([])
return post_summary