使用PyMC的贝叶斯主成分分析

Question

我正在尝试使用Python的PyMC库实现贝叶斯主成分分析（PCA）。但是，我在定义低维坐标时遇到了困难……

模型是

x = Wz + e

这里的x是观察向量，W是转换矩阵，z是低维坐标向量。

首先，我为转换矩阵W定义了一个分布。每一列都是从一个正态分布中抽取的（为了简单起见，均值为零，协方差为单位矩阵）

def W_logp(value):
   logLikes = np.array([multivariate_normal.logpdf(value[:,i], mean=np.zeros(dimX), cov=1) for i in range(0, dimZ)])
   return logLikes.sum()

def W_random():
   W = np.zeros([dimX, dimZ])
   for i in range(0, dimZ):
      W[:,i] = multivariate_normal.rvs(mean=np.zeros(dimX), cov=1)
   return W

w0 = np.random.randn(dimX, dimZ)

W = pymc.Stochastic(
   logp = W_logp,
   doc = 'Transformation',
   name = 'W',
   parents = {},
   random = W_random,
   trace = True,
   value = w0,
   dtype = float,
   rseed = 116.,
   observed = False,
   cache_depth = 2,
   plot = False,
   verbose = 0)

接下来，我想为z定义一个分布，这也是一个多元正态分布（均值为零，协方差为单位矩阵）。不过，我需要为每个观察值单独抽取一个z，而W对所有观察值都是相同的。所以，我尝试了

z = pymc.MvNormal('z', np.zeros(dimZ), np.eye(dimZ), size=N)

然而，pymc.MvNormal没有size参数。这导致了一个错误。接下来的步骤是

m = Data.mean(axis=0) + np.dot(W, z)
obs = pymc.MvNormal('Obs', m, C, value=Data, observed=True)

我没有给出上面C的具体说明，因为现在不相关。有没有什么想法可以实现？

谢谢

编辑

在Chris Fonnesbeck的回答后，我将代码改成了如下

numD, dimX = Data.shape
dimZ = 3
mm = Data.mean(axis=0)

tau = pymc.Gamma('tau', alpha=10, beta=2)
tauW = pymc.Gamma('tauW', alpha=20, beta=2, size=dimZ)

@pymc.deterministic(dtype=float)
def C(tau=tau):
    return (tau)*np.eye(dimX)

@pymc.deterministic(dtype=float)
def CW(tau=tauW):
    return np.diag(tau)

W = [pymc.MvNormal('W%i'%i, np.zeros(dimZ), CW) for i in range(dimX)]
z = [pymc.MvNormal('z%i'%i, np.zeros(dimZ), np.eye(dimZ)) for i in range(numD)]
mu = [pymc.Lambda('mu%i'%i, lambda W=W, z=z: mm + np.dot(np.array(W), np.array(z[i]))) for i in range(numD)]

obs = [pymc.MvNormal('Obs%i'%i, mu[i], C, value=Data[i,:], observed=True) for i in range(numD)]

model = pymc.Model([tau, tauW] + obs + W + z)
mcmc = pymc.MCMC(model)

但这次，当运行pymc.MCMC(model)时，它试图分配大量内存（超过8GB），当numD=45和dimX=504时。即使我只尝试numD=1（只创建1个z、mu和obs），情况也是一样。你知道为什么吗？

暂无标签

使用PyMC的贝叶斯主成分分析

编辑

2 个回答

撰写回答