Gibbs 采样器：数组的第 n 个主要次要不是正定的

Question

我正在尝试构建吉布斯抽样算法来估计线性回归的系数和方差。由于假设方差是常数，我尝试对系数使用 MVNormal 共轭解，对方差使用反伽马共轭解，但这结果非常不准确。

现在我正在尝试实现 MVNormal/Inverse Wishart 解决方案，但我的比例矩阵似乎不是正定的。我收到以下错误：

回溯（最后一次通话）：

文件“D:\STA Classes\STA 141C\Final Project\metropolis_hastings.py”, 第 128 行，在 bayes_betas, bayes_sigmas = gibbs_bayesian(Y, X, np.zeros(100), np.identity(100), 0, np.identity(100), 10000)

文件“D:\STA Classes\STA 141C\Final Project\metropolis_hastings.py”, 第 59 行，在 gibbs_bayesian 中 proposal_sigma = sp.stats.invwishart.rvs(df, psi)

文件 “C:\Users\Sam naconda3\lib\site-packages\scipy\stats_multivariate.py”, 第 2769 行，在 rvs 中 L, lower = scipy.linalg.cho_factor(scale, lower=True)

文件 “C:\Users\Sam naconda3\lib\site-packages\scipy\linalg_decomp_cholesky.py”, 第 152 行，在 cho_factor 中 c, lower = _cholesky(a, lower=lower, overwrite_a=overwrite_a, clean=False,

文件 “C:\Users\Sam naconda3\lib\site-packages\scipy\linalg_decomp_cholesky.py”, 第 37 行，在 _cholesky raise LinAlgError("%d-th leading minor of the array is not positive "

LinAlgError：数组的第 13 个主要次要不是正定的

下面是我的代码：

def gibbs_bayesian(y, X, mu0, sig0, df0, psi0, niter):
    n = X.shape[0]
    p = X.shape[1]
    x_bar = np.sum(X, 0) / p
    beta_estimates = np.empty((p, niter))
    beta_estimates[:,0] = np.ones(p)
    sigma_estimates = np.empty((p, p, niter))
    sigma_estimates[:,:,0] = np.identity(p)
    
    for i in tqdm(range(1, niter)):
        sigma_current = sigma_estimates[:,:,i-1]
        
        # Obtain new estimates of betas
        sig = np.linalg.inv(np.identity(p) + n*np.linalg.inv(sigma_current))
        mu = sig@(np.identity(p)@mu0 + n*np.linalg.inv(sigma_current)@x_bar)
        proposal_beta = np.random.multivariate_normal(mu, sig)
        
        # Obtain new estimates of sigma^2
        df = df0 + n
        rowsums = np.empty((p, p), dtype = "float64")
        for j in range(0, n):
            rowsums += (X[j,:] - proposal_beta) @ (X[j,:] - proposal_beta).T
        psi = psi0 + rowsums
        
        proposal_sigma = sp.stats.invwishart.rvs(df, psi)
        
        # Set new value
        beta_estimates[:,i] = proposal_beta
        sigma_estimates[:,:,i] = proposal_sigma
            
    return beta_estimates, sigma_estimates

我尝试的一件事是向 Psi 矩阵的每个条目添加一个小的正态随机变量。虽然没有运气。任何想法我做错了什么？或者我可以通过另一种方式实现它，从而提供更准确的解决方案？