Created
December 15, 2017 09:59
-
-
Save cdcsai/5c7b3057010acad9b1f2434f9295197b to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
#Estimation, mean of Gamma(4.3, 6.2) with Accept-Reject Algorithm | |
from scipy.stats import gamma | |
import numpy as np | |
import matplotlib.pyplot as plt | |
def f(x): | |
r = gamma.pdf(x, 4.3, scale=1/6.2) | |
return r | |
def accept_reject(n): | |
L = list() | |
u = 3 | |
M = gamma.pdf(u, 4.3, scale=1/6.2) / gamma.pdf(u, 4, scale=1/7) | |
for k in range(n): | |
y = gamma.rvs(4, scale=1/7) | |
u = np.random.uniform() | |
if f(y)/(M*gamma.pdf(y, a=4, scale=1/7)) >= u: | |
x = y | |
L.append(x) | |
return L | |
#Plot Histogram and Density | |
x = np.linspace (-5, 5, 200) | |
L = accept_reject(10000) | |
fig, ax = plt.subplots(1, 1) | |
ax.plot(x, gamma.pdf(x, 4.3, scale=1/6.2), 'r-', label='gamma pdf') | |
ax.hist(L, normed=True) | |
plt.show() | |
#Estimation, mean of Gamma(4.3, 6.2) with Metropolis-Hastings and Gamma(4,6) candidate | |
def target(x): | |
return gamma.pdf(x, 4.3, scale=1 / 6.2) | |
def candidate(x): | |
return gamma.pdf(x, 4, scale=1 / 6) | |
def rho(x, y): | |
return (target(y) / target(x)) * (candidate(x) / candidate(y)) | |
def metropolis_hastings_1(n): | |
samples = np.empty(n) | |
accepted = 0 | |
samples[0] = gamma.rvs(4, scale=1 / 6) | |
for i in range(1, n): | |
xt = samples[i-1] | |
y = gamma.rvs(4, scale=1 / 6) | |
if np.random.uniform() < min(1, rho(xt, y)): | |
samples[i] = y | |
accepted += 1 | |
else: | |
samples[i] = xt | |
return samples | |
x = np.linspace (-5, 5, 200) | |
L = metropolis_hastings_1(10000) | |
fig, ax = plt.subplots(1, 1) | |
ax.plot(x, gamma.pdf(x, 4.3, scale=1/6.2), 'r-', label='gamma pdf') | |
ax.hist(L, normed=True) | |
plt.show() | |
#Metropolis-Hastings, target standard normal distribution and candidate uniform(-x-d, -x+d) | |
from scipy.stats import gamma | |
import numpy as np | |
import matplotlib.pyplot as plt | |
from scipy.stats import norm | |
from scipy.stats import uniform | |
def target(x): | |
return norm.pdf(x, 0, 1) | |
def candidate(x, y, d): | |
return uniform.pdf(y, loc=-x-d, scale=2*d) | |
def rho(x, y, d): | |
return (target(y) / target(x)) * (candidate(y, x, d) / candidate(x, y, d)) | |
def metropolis_hastings_3(n, d): | |
samples = np.empty(n) | |
accepted = 0 | |
x0 = np.random.uniform() | |
samples[0] = uniform.rvs(-x0-d, 2*d) | |
for i in range(1, n): | |
xt = samples[i-1] | |
y = uniform.rvs(-xt-d, 2*d) | |
if np.random.uniform() < min(1, rho(xt, y, d)): | |
samples[i] = y | |
accepted += 1 | |
else: | |
samples[i] = xt | |
return samples | |
x = np.linspace (-5, 5, 200) | |
L = metropolis_hastings_3(10000, 1) | |
fig, ax = plt.subplots(1, 1) | |
ax.plot(x, norm.pdf(x, 0, 1), 'r-', label='gamma pdf') | |
ax.hist(L, normed=True) | |
plt.show() | |
#Gibbs sampler, normal random bivariate variables | |
import numpy as np | |
def gibbs(n, rho): | |
mat = np.zeros((n, 2)) | |
x = 0 | |
y = 0 | |
for i in range(1, n): | |
x = np.random.normal(rho * y, np.sqrt(1 - rho**2)) | |
y = np.random.normal(rho * x, np.sqrt(1 - rho**2)) | |
mat[i][0] = x | |
mat[i][1] = y | |
return mat | |
##Estimation of the X^2+Y^2 density | |
def density_2(n, rho): | |
L = list() | |
R = gibbs(n, rho) | |
for i in range(n): | |
r = (R[i][0])**2 + (R[i][1])**2 | |
L.append(r) | |
return L | |
##Computation of P(X^2+Y^2>2) | |
def indicator(x): | |
if x > 2: | |
x = 1 | |
else: | |
x = 0 | |
return x | |
def proba(n, rho): | |
L = density_2(n, rho) | |
s = 0 | |
for i in range(n): | |
z = np.random.choice(L) | |
s = s + indicator(z) | |
return (1. / n) * s | |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment