自然な共益事前分布を使って線形回帰モデルの係数をベイズ推定

bayesian_liner_regression.py

#! /usr/bin/env python
# -*- coding: utf-8 -*-
"""
自然な共益事前分布を使って線形回帰モデルの係数をベイズ推定
【事前分布】
 平均：多変量正規分布
"""
from numpy.random import *
import numpy as np
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
import sys

def sampling(num_sample,coef1,coef2):
    x1=randn(num_sample)
    x2=randn(num_sample)
    y = coef1 * x1 + coef2 * x2 + randn(num_sample) * 3
    return (x1, x2, y)
    
def estimate_coef(sample_x,sample_y, coef_prior_mean, coef_prior_sigma):
    xx = sample_x.T * sample_x
    coef_mse = xx.I * sample_x.T * sample_y
    coef_sigma = xx + coef_prior_sigma
    coef_mean = (xx + coef_prior_sigma).I * (xx * coef_mse + coef_prior_sigma * coef_prior_mean)
    return (coef_mse,coef_mean)

true_coef1 = 2.
true_coef2 = 4.
prior_coef_mean = np.matrix([2.,4.]).T
prior_coef_sigma = np.matrix(np.eye(2) * 0.1).I
n_sample = 4

np.random.seed(123)
x1, x2, y = sampling(n_sample, true_coef1, true_coef2)
mat_x = np.matrix(zip(x1,x2))

coef_mse, coef_mean = estimate_coef(mat_x,np.matrix(y).T,prior_coef_mean,prior_coef_sigma)
print "最小二乗法:\n{}".format(coef_mse)
print "ベイズ推定:\n{}".format(coef_mean)