Я пытаюсь использовать pyMC для предоставления байесовской оценки ковариационной матрицы с учетом некоторых данных. Я примерно следую примеру фондовой ковариации, приведенному в этом онлайн-руководстве (ссылка здесь), но у меня есть более упрощенная примерная модель, которую я составил. У меня есть два значения, которые я извлек из многомерного нормального распределения, и я построил его таким образом, что знаю ковариацию/корреляцию между двумя переменными.
Я разместил свой короткий код ниже. По сути, я создаю искусственный набор данных, где корреляционная матрица должна быть [[1, -0,5], [-0,5, 1]]. В конце выборки mcmc я получаю прогнозируемое значение недиагонального члена, которое немного отличается. Я посмотрел на критерии сходимости, и похоже, что автокорреляция низкая, а распределение стационарное. Тем не менее, я признаю, что все еще обдумываю все нюансы здесь, и могут быть аспекты этого, которые все еще находятся за пределами моего понимания.
Этот вопрос связан и во многом основан на этих двух других вопросах SO (One и Два). Я чувствовал необходимость задать свой собственный вопрос, несмотря на сходство, потому что я не получаю ответа, который ожидаю получить. Если кто-то из вас, специалистов по вычислительной статистике, может помочь разобраться в этой проблеме, мы будем очень признательны!
import numpy as np
import pandas as pd
import pymc as pm
import matplotlib.pyplot as plt
import seaborn as sns
p=2
prior_mu=np.ones(p)
prior_sdev=np.ones(p)
prior_corr_inv=np.eye(p)
def cov2corr(A):
"""
covariance matrix to correlation matrix.
"""
d = np.sqrt(A.diagonal())
A = ((A.T / d).T) / d
#A[ np.diag_indices(A.shape[0]) ] = np.ones( A.shape[0] )
return A
# construct artificial data set
muVector=[10,5]
sdevVector=[3.,5.]
corrMatrix=np.matrix([[1,-0.5],[-0.5, 1]])
cov_matrix=np.diag(sdevVector)*corrMatrix*np.diag(sdevVector)
n_obs = 500
x = np.random.multivariate_normal(muVector,cov_matrix,n_obs)
prior_mu = np.array(muVector)
prior_std = np.array(sdevVector)
inv_cov_matrix = pm.Wishart( "inv_cov_matrix", n_obs, np.diag(prior_std**2) )
mu = pm.Normal( "returns", prior_mu, 1, size = 2)
# create the model and sample
obs = pm.MvNormal( "observed returns", mu, inv_cov_matrix, observed = True, value = x )
model = pm.Model( [obs, mu, inv_cov_matrix] )
mcmc = pm.MCMC(model)
mcmc.use_step_method(pm.AdaptiveMetropolis,inv_cov_matrix)
mcmc.sample( 1e5, 2e4, 10)
# Determine prediction - Does not equal corrMatrix!
inv_cov_samples = mcmc.trace("inv_cov_matrix")[:]
mean_covariance_matrix = np.linalg.inv( inv_cov_samples.mean(axis=0) )
prediction = cov2corr(mean_covariance_matrix*n_obs)
