在pythone中从R重新创建Gibbs采样器

Question

下面是我收到的错误。我对Python编程相当陌生。任务是用不同于R的语言重新创建Gibbs采样器，因此我选择Python来帮助自己熟悉该语言。非常感谢您对这个错误的任何帮助

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-156-56d5c35de2b3> in <module>()
 50 x2_loc = data.drop(columns = ['y', 'Unnamed: 0'])
 51 
---> 52 gibbs_horseshoe(y_loc, x2_loc, 100)

<ipython-input-156-56d5c35de2b3> in gibbs_horseshoe(y, X, N_sims)
 33     #begin for loop
 34     for i in range(1, N_sims + 1):
---> 35         SIGMA = np.linalg.inv(XtX + np.zeros((tau_loc, p)))
 36         beta_loc = beta_out[:, i + 1] = 
np.array([np.random.normal(np.dot(SIGMA, Xty), sigma_sq_loc*SIGMA, 1)])
 37         sigma_sq_loc = sigma_sq_out[i + 1] = 1/np.random.gamma((n + p + 
1)/2, (1/2)*np.dot(np.dot(y-X, beta_loc)) + np.dot(tau_loc*beta_loc, 
beta_loc) + 1, 1)

TypeError: only integer scalar arrays can be converted to a scalar index
---------------------------------------------------------------------------

代码如下

import numpy as np
import pandas as pd

def gibbs_horseshoe(y, X, N_sims):


Xt = np.transpose(X)
XtX = np.dot(Xt,X)
Xty = np.dot(Xt,y)
p = np.shape(X)[1] #Number of columns
n = np.shape(X)[0] #Number of rows
N_sims = 1

tau_out = beta_out = np.zeros((p, N_sims + 1))
sigma_sq_out = np.zeros((N_sims + 1))

tau_loc = lambda_loc = np.ones(p)
phi_loc = eta_loc = 1

diag_p = np.zeros((p, p))

beta_loc = np.dot(np.linalg.inv(XtX + np.fill_diagonal(diag_p, p)/n), Xty)
sigma_sq_loc = 1/n(np.cross(y, y) - np.cross(Xty, beta_loc))[1, 1]

tau_out[:, 1] = tau_loc
beta_out[:, 1] = beta_loc
sigma_sq_out[1] = sigma_sq_loc

#begin for loop
for i in range(1, N_sims + 1):
    SIGMA = linalg.inv(XtX + np.zeros(tau_loc, p))
    beta_loc = beta_out[:, i + 1] = np.array([np.random.normal(np.dot(SIGMA, Xty), sigma_sq_loc*SIGMA, 1)])
    sigma_sq_loc = sigma_sq_out[i + 1] = 1/np.random.gamma((n + p + 1)/2, (1/2)*np.cross(np.dot(y-X, beta_loc)) + np.cross(tau_loc*beta_loc, beta_loc) + 1, 1)
    tau_loc = tau_out[:, i + 1] = np.random.exponential(1, p)/(beta_loc^2/(2*sigma_sq_loc)+lambda_loc/2)
    lambda_loc = np.random.exponential(1, p)/(tau_loc/2+phi_loc/2)
    phi_loc = np.random.gamma((p+1)/2, sum(lambda_loc)/2+eta_loc/2, 1)
    eta_loc = np.random.exponential(phi_loc/2+1/2, 1)
#end for loop


return (beta_out, sigma_sq_out, tau_out)

filename = 'diabetes3.csv'
data = pd.read_csv(filename)
y_loc = data['y']
x2_loc = data.drop(columns = ['y', 'Unnamed: 0'])

gibbs_horseshoe(y_loc, x2_loc, 100)

Answer 1

np.zeros(tau_loc, p) 这部分既不起作用也没有意义

我真的不认为我应该查找所有的吉布斯采样代码。 但我认为你应该从那里开始工作