PyMC - 方差 - 協方差矩陣估計

我閱讀以下紙（http://www3.stat.sinica.edu.tw/statistica/oldpdf/A10n416.pdf），其中它們的方差 - 協方差矩陣Σ建模爲：PyMC - 方差 - 協方差矩陣估計

Σ= DIAG（S）* R * DIAG（S）...（式1中（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)

來源

2014-02-10 akhil

請詳細解釋「模型」的含義。這個詞在統計和科學中有很多含義，這些在這裏似乎都不適用。你可能在問怎麼把協方差矩陣分解成這種形式？如果你的問題只是關於在PyMC中編碼算法，那麼請讓我們知道，以便我們可以將它遷移到SO社區。 – whuber

我的問題只是關於PyMC中的實現。 – akhil

我認爲這可以使用以下方法：@ 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） – akhil

這裏該解決方案有益於某人絆倒這個職位：

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)

來源

2014-02-11 19:26:57 akhil

PyMC - 方差 - 協方差矩陣估計

回答

相關問題