自然な共益事前分布を使って正規分布のパラメーターを推定。平均：正規分布、分散:逆ガンマ分布。

bayesian_normal_dist.py

#! /usr/bin/env python
# -*- coding: utf-8 -*-
"""
自然な共益事前分布を使って正規分布のパラメーターをベイズ推定
【事前分布】
 平均：正規分布
 分散:逆ガンマ分布
"""
import numpy as np
import scipy.stats
import scipy.special
from numpy.random import *
import matplotlib.pyplot as plt
import math
import sys



def sampling(num_sample,mean=0.,scale=1.):
    sample = randn(num_sample) * np.sqrt(scale) + mean
    return sample

def estimate_mean(sample,prior_mean,prior_sigma,true_sigma):
    """
    事前分布に正規分布を用いて平均を推定する
    """
    n = float(sample.shape[0])
    posterior_sigma = 1./ (1. / prior_sigma + n / prior_sigma)
    posterior_mean = sigma / (n * prior_sigma + true_sigma) * prior_mean + prior_sigma / (n * prior_sigma + true_sigma) * sample.sum()
    return (posterior_mean,posterior_sigma)
    

def invGamma(x,alpha,beta):
    return (beta ** alpha / scipy.special.gamma(alpha)) * x ** - (alpha + 1) * np.exp(-beta/x)
    
def estimate_sigma(sample,prior_alpha,prior_beta,true_mean):
    """
    事前分布に逆ガンマ分布を用いて分散を推定する
    """
    n = float(sample.shape[0])
    sample_var = n * np.var(sample)
    alpha = (n + prior_alpha) / 2
    beta = (prior_beta + sample_var) / 2
    
    return (beta / (alpha - 1), beta ** 2 / ((alpha - 1) ** 2 * (alpha - 2)))

def estimate_mean_var(sample, mean_prior_mean, mean_prior_sigma, sigma_prior_alpha, sigma_prior_beta):
    """
    事前分布に正規分布と逆ガンマ分布を用いて平均と分散を推定する
    """
    n = float(sample.shape[0])
    alpha = (sigma_prior_alpha + n) / 2
    beta = (sigma_prior_beta + (n - 1) * np.var(sample) + (mean_prior_sigma * n) / (mean_prior_sigma + n) * (mean_prior_mean - np.average(sample))) / 2
    sigma_posterior_mean = beta / (alpha - 1)
    sigma_posterior_sigma = beta ** 2 / ((alpha - 1) ** 2 * (alpha - 2))
    mean_posterior_mean = mean_prior_sigma / (mean_prior_sigma + n) * mean_prior_mean + n / (mean_prior_sigma + n) * np.average(sample)
    mean_posterior_sigma = sigma_posterior_mean / (mean_prior_sigma + n)
    return (mean_posterior_mean,mean_posterior_sigma,sigma_posterior_mean,sigma_posterior_sigma)

if __name__ == "__main__":
    #サンプルの真のパラメーター(平均と分散)
    mean = 20.
    sigma = 11. 
    num_sample = 4

    np.random.seed(100)

    sample_train = sampling(num_sample,mean=mean,scale=sigma)

    ############# 平均を推定 #############
    #事前分布
    mean_prior_mean = 20.
    mean_prior_sigma = 0.1

    posterior_mean, posterior_sigma = estimate_mean(sample_train,mean_prior_mean,mean_prior_sigma,sigma)

    print u"平均の最尤推定値:%f"%(np.average(sample_train))
    print u"平均のベイズ推定値:%f"%(posterior_mean)

    x = np.linspace(18,22,100)
    y = scipy.stats.norm.pdf(x,loc=posterior_mean,scale=np.sqrt(posterior_sigma))
    plt.title('posterior distribution of mean')
    plt.plot(x,y)
    plt.show()

    ############# 分散を推定 #############

    sigma_prior_mean = 10.
    sigma_prior_sigma = 1.
    # Reference: http://blog.livedoor.jp/incaseforgetr/archives/40094996.html
    #　check:２で割るか否か
    prior_alpha=sigma_prior_mean ** 2 / sigma_prior_sigma-2
    prior_beta=np.sqrt( sigma_prior_sigma *( sigma_prior_mean ** 2 / sigma_prior_sigma -3  ) ** 2    *   (sigma_prior_mean ** 2 /sigma_prior_sigma -4))
    posterior_mean, _ = estimate_sigma(sample_train,prior_alpha,prior_beta,mean)
    print u"分散の最尤推定値:%f"%(np.var(sample_train))
    print u"分散のベイス推定値:%f"%(posterior_mean)

    ############# 平均と分散を推定 #############
    mean_mean, _ , sigma_mean, _ = estimate_mean_var(sample_train,mean_prior_mean,mean_prior_sigma,prior_alpha ,prior_beta)
    print u"平均と分散のベイズ推定:{}".format((mean_mean,sigma_mean))