PyMC - 方差-协方差矩阵估计

Posted 2023-02-19

【中文标题】PyMC - 方差-协方差矩阵估计

【英文标题】：PyMC - variance-covariance matrix estimation

【发布时间】：2014-03-09 18:54:33

【问题描述】：

我阅读了以下论文(http://www3.stat.sinica.edu.tw/statistica/oldpdf/A10n416.pdf)，他们将方差-协方差矩阵 Σ 建模为：

Σ = diag(S)*R*diag(S)（论文中的公式1）

S是k×1的标准差向量，diag(S)是对角元素S的对角矩阵，R是k×k相关矩阵。

如何使用 PyMC 来实现？

这是我写的一些初始代码：

import numpy as np
import pandas as pd
import pymc as pm

k=3
prior_mu=np.ones(k)
prior_var=np.eye(k)
prior_corr=np.eye(k)
prior_cov=prior_var*prior_corr*prior_var

post_mu = pm.Normal("returns",prior_mu,1,size=k)
post_var=pm.Lognormal("variance",np.diag(prior_var),1,size=k)
post_corr_inv=pm.Wishart("inv_corr",n_obs,np.linalg.inv(prior_corr))


post_cov_matrix_inv = ???

muVector=[10,5,-2]
varMatrix=np.diag([10,20,10])
corrMatrix=np.matrix([[1,.2,0],[.2,1,0],[0,0,1]])
cov_matrix=varMatrix*corrMatrix*varMatrix

n_obs=10000
x=np.random.multivariate_normal(muVector,cov_matrix,n_obs)
obs = pm.MvNormal( "observed returns", post_mu, post_cov_matrix_inv, observed = True, value = x )

model = pm.Model( [obs, post_mu, post_cov_matrix_inv] )
mcmc = pm.MCMC()

mcmc.sample( 5000, 2000, 3 )

谢谢

[编辑]

我认为可以使用以下方法来完成：

@pm.deterministic
def post_cov_matrix_inv(post_sdev=post_sdev,post_corr_inv=post_corr_inv):
    return np.diag(post_sdev)*post_corr_inv*np.diag(post_sdev)

【问题讨论】：

请详细说明“模型”的含义。这个词在统计学和科学中有很多含义，但似乎都不适用于这里。您是否在问如何将协方差矩阵分解成这种形式？如果您的问题只是关于在 PyMC 中编写算法，请告诉我们，以便我们将其迁移到 SO 社区。​​span>

我的问题只是关于 PyMC 中的实现。

我认为可以使用以下方法来完成：@pm.deterministic def post_cov_matrix_inv(post_sdev=post_sdev,post_corr_inv=post_corr_inv): return np.diag(post_sdev)*post_corr_inv*np.diag(post_sdev)

【参考方案1】：

以下是偶然发现此帖子的人的解决方案：

p=3
prior_mu=np.ones(p)
prior_sdev=np.ones(p)
prior_corr_inv=np.eye(p)


muVector=[10,5,1]
sdevVector=[3,5,10]
corrMatrix=np.matrix([[1,0,-.1],[0,1,.5],[-.1,.5,1]])
cov_matrix=np.diag(sdevVector)*corrMatrix*np.diag(sdevVector)

n_obs=2000
x=np.random.multivariate_normal(muVector,cov_matrix,n_obs)

prior_cov=np.diag(prior_sdev)*np.linalg.inv(prior_corr_inv)*np.diag(prior_sdev)

post_mu = pm.Normal("returns",prior_mu,1,size=p)
post_sdev=pm.Lognormal("sdev",prior_sdev,1,size=p)
post_corr_inv=pm.Wishart("inv_corr",n_obs,prior_corr_inv)

#post_cov_matrix_inv = pm.Wishart("inv_cov_matrix",n_obs,np.linalg.inv(prior_cov))
@pm.deterministic
def post_cov_matrix_inv(post_sdev=post_sdev,post_corr_inv=post_corr_inv,nobs=n_obs):
    post_sdev_inv=(post_sdev)**-1
    return np.diag(post_sdev_inv)*cov2corr(post_corr_inv/nobs)*np.diag(post_sdev_inv)

obs = pm.MvNormal( "observed returns", post_mu, post_cov_matrix_inv, observed = True, value = x )

model = pm.Model( [obs, post_mu, post_sdev ,post_corr_inv])
mcmc = pm.MCMC(model)

mcmc.sample( 25000, 15000, 1,progress_bar=False )