various computational stats algorithms.py

#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