<jupyter_start><jupyter_text>Gaussian inferences<jupyter_code>data = np.array([51.06, 55.12, 53.73, 50.24, 52.05, 56.40, 48.45, 52.34, 55.65, 51.49, 51.86, 63.43, 53.00, 56.09, 51.93, 52.31, 52.33, 57.48, 57.44, 55.14, 53.93, 54.62, 56.09, 68.58, 51.36, 55.47, 50.73, 51.94, 54.95, 50.39, 52.91, 51.5, 52.68, 47.72, 49.73, 51.82, 54.99, 52.84, 53.19, 54.52, 51.46, 53.73, 51.61, 49.81, 52.42, 54.3, 53.84, 53.16])

# remove outliers using the interquartile rule
quant = np.percentile(data, [25, 75])
iqr = quant[1] - quant[0]
upper_b = quant[1] + iqr * 1.5
lower_b = quant[0] - iqr * 1.5
clean_data = data[(data > lower_b) & (data < upper_b)]
sns.kdeplot(data)
plt.xlabel('$x$', fontsize=16)
with pm.Model() as model_g:
    mu = pm.Uniform('mu', lower=40, upper=70)
    sigma = pm.HalfNormal('sigma', sd=10)
    y = pm.Normal('y', mu=mu, sd=sigma, observed=data)
    trace_g = pm.sample(1100, njobs=1)<jupyter_output>Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Sequential sampling (2 chains in 1 job)
NUTS: [sigma, mu]
100%|██████████| 1600/1600 [00:01<00:00, 1357.11it/s]
100%|██████████| 1600/1600 [00:00<00:00, 2131.76it/s]<jupyter_text>Note that the traces of each variable are the marginal distributions not the real posterior which is a 2-dimensionsl joint distribution<jupyter_code>chain_g = trace_g[100:]
pm.traceplot(chain_g)
pm.summary(chain_g)
y_pred = pm.sample_ppc(chain_g, 100, model_g, size=len(data))
sns.kdeplot(data, color='b')
for i in y_pred['y']:
    sns.kdeplot(i.flatten(), color='r', alpha=0.1)
plt.title('Gaussian model', fontsize=16)
plt.xlabel('$x$', fontsize=16)<jupyter_output>/home/damianos/miniconda3/envs/pymc3/lib/python3.7/site-packages/ipykernel_launcher.py:1: DeprecationWarning: sample_ppc() is deprecated.  Please use sample_posterior_predictive()
  """Entry point for launching an IPython kernel.
100%|██████████| 100/100 [00:00<00:00, 1442.89it/s]<jupyter_text>Gaussian Robust inferencesNow we are going to learn how to make a robust model to outliers, that is a model that do not get too much excited by outliers. For that we are going to use a Student T distribution which can have heavier tails and allows for the presence of outliers<jupyter_code>plt.figure(figsize=(8, 6))
x_values = np.linspace(-10, 10, 200)
for df in [1, 2, 5, 30]:
    distri = stats.t(df)
    x_pdf = distri.pdf(x_values)
    plt.plot(x_values, x_pdf, label=r'$\nu$ = {}'.format(df))

x_pdf = stats.norm.pdf(x_values)
plt.plot(x_values, x_pdf, label=r'$\nu = \infty$')
plt.xlabel('$x$', fontsize=16)
plt.ylabel('$pdf(x)$', fontsize=16, rotation=90)
plt.legend(loc=0, fontsize=14)
plt.xlim(-7, 7);
with pm.Model() as model_t:
    mu = pm.Uniform('mu', 40, 75)
    sigma = pm.HalfNormal('sigma', sd=10)
    nu = pm.Exponential('nu', 1/30)
    y = pm.StudentT('y', mu=mu, sd=sigma, nu=nu, observed=data)
    trace_t = pm.sample(11000,njobs=8)
chain_t = trace_t[100:]
pm.traceplot(chain_t);
pm.summary(chain_t)<jupyter_output><empty_output><jupyter_text>we can see that $\nu$~4 which indicates heavier tails than the normal when $\nu$->inf<jupyter_code>y_pred = pm.sample_ppc(chain_t, 100, model_t, size=len(data))
sns.kdeplot(data, c='b')
for i in y_pred['y']:
    sns.kdeplot(i.flatten(), c='r', alpha=0.1)
plt.xlim(35, 75)
plt.title("Student's t model", fontsize=16)
plt.xlabel('$x$', fontsize=16)<jupyter_output>/home/damianos/miniconda3/envs/pymc3/lib/python3.7/site-packages/ipykernel_launcher.py:1: DeprecationWarning: sample_ppc() is deprecated.  Please use sample_posterior_predictive()
  """Entry point for launching an IPython kernel.
100%|██████████| 100/100 [00:00<00:00, 1081.86it/s]
/home/damianos/miniconda3/envs/pymc3/lib/python3.7/site-packages/matplotlib/cbook/__init__.py:1725: UserWarning: Saw kwargs ['c', 'color'] which are all aliases for 'color'.  Kept value from 'color'
  seen=seen, canon=canonical, used=seen[-1]))<jupyter_text>Tips example  we want to investigate the effect the day of the week has on the tips<jupyter_code>tips = sns.load_dataset('tips')
tips.tail()
sns.violinplot(x='day', y='tip', data=tips);
y = tips['tip'].values
x = pd.Categorical(tips['day']).codes
categories = pd.Categorical(tips['day']).categories
print('Tips sample: {}'.format(y[:5]))
print('Encoded days sample: {}'.format(x[:5]))
print('Decoded days sample: {}'.format(categories[x[:5]]))
with pm.Model() as comparing_groups:    
    # priors
    means = pm.Normal('means', mu=0, sd=10, shape=len(set(x)))
    sds = pm.HalfNormal('sds', sd=10, shape=len(set(x)))
    # a different trace for each weekday as encoded in `x`
    y_pred = pm.Normal('y', mu=means[x], sd=sds[x], observed=y)    
    trace_cg = pm.sample(5000, chains=8, njobs=8)
chain_cg = trace_cg[100::]
pm.traceplot(chain_cg);
pm.summary(chain_cg)
dist = dist = stats.norm()

_, ax = plt.subplots(3, 2, figsize=(16, 12))

comparisons = [(i,j) for i in range(4) for j in range(i+1, 4)]
pos = [(k,l) for k in range(3) for l in (0, 1)]

print('Comparison between the followin encoded days: {}'.format(comparisons))
print('Subplot indices: {}'.format(pos))
for (i, j), (k,l) in zip(comparisons, pos):
    # mean difference between days to be compared
    means_diff = chain_cg['means'][:,i]-chain_cg['means'][:,j]
    # effect size via Cohen's delta
    d_cohen = (means_diff / np.sqrt((chain_cg['sds'][:,i]**2 + chain_cg['sds'][:,j]**2) / 2)).mean()
    # probability of superiority
    ps = dist.cdf(d_cohen/(2**0.5))
        
    pm.plot_posterior(means_diff, ref_val=0, ax=ax[k, l], kde_plot=True, lw=2)
    ax[k, l].plot(0, label="Cohen's d = {:.2f}\nProb sup = {:.2f}".format(d_cohen, ps) ,alpha=0)
    ax[k, l].set_xlabel('$\mu_{}-\mu_{}$ ({} - {})'.format(i, j,categories[i],categories[j]), fontsize=18)
    ax[k,l ].legend(loc=0, fontsize=14)

    
plt.tight_layout()<jupyter_output>Comparison between the followin encoded days: [(0, 1), (0, 2), (0, 3), (1, 2), (1, 3), (2, 3)]
Subplot indices: [(0, 0), (0, 1), (1, 0), (1, 1), (2, 0), (2, 1)]<jupyter_text>We can see that the only comparison where HPD95% excludes the reference values of zero is Thursday and Sunday. For the rest we cannot rule out a difference of zero. But the mean difference for Thursday - Sundsay is ~0.5, is this worth working weekends? Statisticians cannot answer these kind of questions For Cohen's delta see https://en.wikipedia.org/wiki/Effect_sizeCohen's_d  Probability of superiority $ps=\Phi(\frac{\delta}{\sqrt(2)})$ if we assume the group disrtibutions are Normal (not really used much)  Hierarchical Models<jupyter_code>N_samples =  [30, 30, 30]
# Number of good samples, i.e samples with lead content below WHO recommendations
G_samples =  [18, 3, 3] # [13, 3, 3] [18, 3, 3] [18, 18, 18]

group_idx = np.repeat(np.arange(len(N_samples)), N_samples)
data = []  
for i in range(0, len(N_samples)):
    data.extend(np.repeat([1, 0], [G_samples[i], N_samples[i]-G_samples[i]]))
print('Group id of each sample: {}'.format(group_idx))
print('Samples : (0->Bad, 1->Good) {}'.format(data))
print('Total number of samples = {}'.format(len(data)))<jupyter_output>Group id of each sample: [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1
 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2
 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2]
Samples : (0->Bad, 1->Good) [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
Total number of samples = 90<jupyter_text>$$\alpha \sim HalfCauchy(\beta_{\alpha})$$$$\beta \sim HalfCauchy(\beta_{\beta})$$$$\theta \sim Beta(\alpha, \beta)$$$$y \sim Bern(\theta)$$<jupyter_code>with pm.Model() as model_h:
    # hyper-priors
    alpha = pm.HalfCauchy('alpha', beta=10)
    beta = pm.HalfCauchy('beta', beta=10)
    # prior
    theta = pm.Beta('theta', alpha, beta, shape=len(N_samples))
    # likelihood
    y = pm.Bernoulli('y', p=theta[group_idx], observed=data)

    trace_h = pm.sample(8000,njobs=8,chains=8)
chain_h = trace_h[200:]
pm.traceplot(chain_h)
pm.summary(chain_h)
x = np.linspace(0, 1, 100)
print(chain_h['theta'].shape)

for i in np.random.randint(0, len(chain_h), size=100):
    pdf = stats.beta(chain_h['alpha'][i], chain_h['beta'][i]).pdf(x)
    plt.plot(x, pdf,  'g', alpha=0.1)

dist = stats.beta(chain_h['alpha'].mean(), chain_h['beta'].mean())
pdf = dist.pdf(x)
mode = x[np.argmax(pdf)]
mean = dist.moment(1)
plt.plot(x, pdf, label='mode = {:.2f}\nmean = {:.2f}'.format(mode, mean))
plt.title(r'Dispersion of the distribution of $\theta$ based on the $\alpha, \beta$ hyperpriors')

plt.legend(fontsize=14)
plt.xlabel('$\\theta_{prior}$', fontsize=16)
plt.tight_layout()

plt.figure()
for i, thetaChain in enumerate(chain_h['theta'].T):
    plt.hist(thetaChain,alpha=0.3, label='Probability of water being good for group {}'.format(i))
plt.legend()
import sys, IPython, scipy, matplotlib, platform
print("This notebook was created on a %s computer running %s and using:\nPython %s\nIPython %s\nPyMC3 %s\nNumPy %s\nSciPy %s\nMatplotlib %s\nSeaborn %s\nPandas %s" % (platform.machine(), ' '.join(platform.linux_distribution()[:2]), sys.version[:5], IPython.__version__, pm.__version__, np.__version__, scipy.__version__, matplotlib.__version__, sns.__version__, pd.__version__))<jupyter_output>This notebook was created on a x86_64 computer running debian buster/sid and using:
Python 3.7.2
IPython 7.2.0
PyMC3 3.6
NumPy 1.16.0
SciPy 1.2.0
Matplotlib 3.0.2
Seaborn 0.9.0
Pandas 0.23.4
{'fraction_non_alphanumeric': 0.1826712543009071, 'fraction_numerical': 0.07882389740381608, 'mean_word_length': 4.629107981220657, 'pattern_counts': {'":': 0, '<': 28, '<?xml version=': 0, '>': 29, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 19, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '185639', 'license': 'MIT', 'n_tokens_mistral': 4498, 'n_tokens_neox': 3844, 'n_words': 1054, 'path': 'Chapter 3/Bayesian Data Analysis Chapter 3.ipynb', 'repo_name': 'chdamianos/Bayesian-Analysis-with-Python', 'size': '9591'}