CoCalc Public Filesatms391geodata / Week 11 / Homework 11 (answer key).ipynb
Author: Steve Nesbitt
Views : 51
Compute Environment: Ubuntu 18.04 (Deprecated)

ATMS 391 Homework 11

Problem 1

Assume that you have a hypothetical situation in which it is claimed that the climatological probability of a cloudless day in winter is 6/7, after observing x = 15 cloudless days on N = 25 independent occasions. Assume that the prior distribution can be approximated by a Beta distrubution with $\alpha = \beta =$ 4.

(a) Generate the Posterior distribution using Bayes' theorem.

In [13]:
import numpy as np
import matplotlib.pyplot as plt
from __future__ import division
%matplotlib inline
import scipy.stats as st

def norm_beta(n_trial, n_cloud, alpha_prior, beta_prior):
x_values = np.linspace(0, 1, 100)
prior = [np.product(st.beta.pdf(x1, alpha_prior, beta_prior)) for x1 in x_values]
p = 15 / 25 # probability
mu = n_trial * p
sigma = (n_trial * p * (1 - p))**0.5
likelihood = [np.product(st.norm.pdf(x1, loc=mu, scale=sigma)) for x1 in x_values]
posterior = [st.beta.pdf(x, alpha_prior + n_cloud, beta_prior + (n_trial - n_cloud)) for x in x_values]

fig, axes = plt.subplots(3, 1, sharex=True, figsize=(8,8))
axes[0].plot(x_values, likelihood)
axes[0].set_title("Sampling Distribution")
axes[1].plot(x_values, prior)
axes[1].set_title("Prior Distribution")
axes[2].plot(x_values, posterior)
axes[2].set_title("Posterior Distribution")
plt.tight_layout()
return posterior, x_values

posterior, x_values = norm_beta(25, 15, 4, 4)


(b) What is the mean of the posterior distribution?

In [14]:
n = 25
y = 15
prior_alpha = 4
prior_beta = 4
posterior_mean = (y + prior_alpha) / (y + prior_alpha + n - y + prior_beta)
print("Posterior Mean: {}".format(posterior_mean))

Posterior Mean: 0.575757575758

(c) What are the 95% confidence intervals of the Posterior?

In [23]:
conf_int = np.percentile(posterior, [2.5,97.5])
conf_int

[7.5531402197569343e-22, 4.5411267644511275]

(d) What is the probability that the original claim is correct? (Hint: Evaluate the sum of the posterior above the probability 6/7)

In [19]:
pdf = posterior/np.sum(posterior)

probability= np.sum(pdf[x_values > 6/7])
probability

6.9691911785028847e-05

Problem 2

Run Exercise 4, we will talk about it next Tuesday!!

In [ ]: