首页 文章

使用PyMC的贝叶斯pca

提问于
浏览
2

我正在尝试使用PyMC库为python实现贝叶斯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) 时使用 numD=45dimX=504 分配大量内存(超过8GB) . 即使我只用 numD=1 (因此仅创建1 zmuobs )尝试它,它也会这样做 . 知道为什么吗?

2 回答

  • 3

    不幸的是,PyMC不容易让你定义多变量随机指标的向量 . 希望我们可以在PyMC 3中实现这一点 . 现在,您必须使用容器指定它 . 例如:

    z = [pymc.MvNormal('z_%i' % i, np.zeros(dimZ), np.eye(dimZ)) for i in range(N)]
    
  • 1

    关于内存问题,请尝试使用不同的后端作为跟踪 . 默认值( "ram" )将所有内容保存在RAM中 . 您可以尝试使用 "pickle""sqlite" 之类的内容 .

    关于板块符号,它可能是我们可以为PyMC 3寻求的东西 . 随意在_2435558中创建一个问题 .

相关问题