CoCalc Public Filesatms391geodata / Week 11 / Homework 11 (answer key).ipynbOpen with one click!
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 [ ]: