notebook.community

Edit and run



In [28]:

    
import numpy as np
from numpy import linalg as la
import pandas as pd

ギブスサンプラーによる回帰モデルの推定

①ベイズ統計の枠組みに基づく回帰モデルの表現

K個の独立変数$$x_k(k=1, ..., K)$$によって従属変数yを説明する回帰分析のモデル式は、以下のように表すことができる。$$y_i=\beta_1x_{1i}+ ... +\beta_Kx_{Ki}+\epsilon_i$$

標準的なベイズ統計による回帰分析では、誤差項εiの分布がN(0, σ^2)に従うものと仮定する。この時、上の式のy_iは、各β_kx_kiにランダムな正規誤差ε_iが加わって生成されるものと見なすことができる。

私は



In [5]:

    
x_1 = np.ones(10)
x_2 = np.zeros(10)
for i in range(10):
    x_2[i] = np.random.normal(0, 2)
y = 3*x_1+(-2)*x_2



In [15]:

    
pd.DataFrame({"x_1":x_1, "x_2": x_2, "y":y})



In [ ]:

    
beta_given_var =



In [ ]:

    
class GibbsSampler:



In [25]:

    
X = np.matrix([x_1, x_2]).T



In [26]:

    
X.T*X









    Out[26]:





matrix([[ 10.        , -10.37791064],
        [-10.37791064,  71.05328658]])



In [32]:

    
X_t_and_X = X.T*X



In [29]:

    
beta = np.matrix([[1000, 0], [0, 1000]])



In [31]:

    
beta_inv = la.inv(beta)



In [35]:

    
beta_inv









    Out[35]:





matrix([[ 0.001,  0.   ],
        [ 0.   ,  0.001]])



In [37]:

    
beta_2_inv = beta_inv + X_t_and_X



In [38]:

    
beta_2 = la.inv(beta_2_inv)



In [39]:

    
beta_2









    Out[39]:





matrix([[ 0.11785166,  0.01721295],
        [ 0.01721295,  0.0165878 ]])



In [ ]:

	x_1	x_2	y
0	1	0.955671	1.088657
1	1	-3.193369	9.386739
2	1	3.056116	-3.112232
3	1	-4.043210	11.086420
4	1	2.006564	-1.013127
5	1	0.550046	1.899908
6	1	-3.823734	10.647468
7	1	-3.535172	10.070345
8	1	-1.026138	5.052276
9	1	-1.324685	5.649369