多元d拟合模型比较

Question

我在PyMC3中使用WAIC（广泛适用的信息准则）时遇到问题。也就是说，我知道数据是按照多元Dirichlet分布来分布的。我试图通过假设边际分布在一种情况下是贝塔分布，而在另一种情况下是对数正态分布来拟合数据。显然，在第一种情况下，我得到的WAIC值比第二种情况低（更好）。在

第三种情况下，根据Dirichlet分布，假设问题出现在第三种情况下。第三个WAIC明显比前两个病例大（差）。我希望这个WAIC比我在第二个（对数正常）情况下得到的更低（更好）。在

基本上我想证明对数正态拟合是不好的。这是肉眼很容易看到的，但我想有正式的结果来显示。在

复制我得到的东西的最小代码：

import pandas as pd
import numpy as np
import pymc3 as pm

# generate the data
df=pd.DataFrame(np.random.dirichlet([10,10,10],size=2000))

# fit the first case (assuming beta marginal distributions)
betaModel=pm.Model()
with betaModel:
    alpha=pm.Uniform("alpha",lower=0,upper=20,shape=3)
    beta=pm.Uniform("beta",lower=0,upper=40,shape=3)
    observed=pm.Beta("obs",alpha=alpha,beta=beta,observed=df.values,shape=df.shape)
    betaTrace=pm.sample()

# fit the second case (assuming log-normal marginal distributions)
lognormalModel=pm.Model()
with lognormalModel:
    mu=pm.Normal("mu",mu=0,sd=3,shape=3)
    sd=pm.HalfNormal("sd",sd=3,shape=3)
    observed=pm.Lognormal("obs",mu=mu,sd=sd,observed=df.values,shape=df.shape)
    lognormalTrace=pm.sample()

# fit the third case (assuming Dirichlet multivariate distribution)
dirichletModel=pm.Model()
with dirichletModel:
    alpha=pm.HalfNormal("alpha",sd=3,shape=3)
    observed=pm.Dirichlet("obs",a=alpha,observed=df.values,shape=df.shape)
    dirichletTrace=pm.sample()

# compare WAIC
print(pm.waic(betaTrace,betaModel))
print(pm.waic(lognormalTrace,lognormalModel))
print(pm.waic(dirichletTrace,dirichletModel))

输出为：

^{2}$

我想问题可能与错误有关：

ValueError: operands could not be broadcast together with shapes (6000,) (2000,) 

当我试着跑的时候

pm.compare((betaTrace,lognormalTrace,dirichletTrace),(betaModel,lognormalModel,dirichletModel))

如何进行合理的比较有什么建议？在

编辑

在考虑了这个问题之后，我认为这有点“不恰当”。我倾向于这样认为，因为WAIC是一个相对的度量，因此很可能只有相似的统计模型才能进行合理的比较。如果模型太不一样，那你就得到我的了。在

我从pm.compare得到的错误似乎与如何处理随机向量有关。在前两种情况下，随机向量的每个分量都被视为单独的随机变量（每2000个向量有3个分量=6000个点）。在第三种情况下，整个随机向量被视为一个随机变量（2000个向量=2000个点）。在

一开始我认为这个问题可以通过减少前两种情况下的点数来解决。但是，由于前两个统计模型（错误地）假设组件是独立的，增加日志概率不会改变任何事情。WAIC值保持不变。在

目前我认为一个小小的欺骗是可能的。也就是说，将数据拟合到Dirichlet分布中，但是计算WAIC就像我会拟合beta分布一样。这给出了一个预期的结果-Dirichlet拟合的WAIC略大于beta拟合的WAIC，但小于对数正态拟合的WAIC。在

这个“作弊”的密码是：

from collections import namedtuple
from scipy.special import logsumexp

def cheat_logp(tracePoint,model):
    values=model.obs.eval()
    _,components=values.shape
    cb=[None]*components
    beta=np.sum(tracePoint["alpha"])
    for i in range(components):
        cheatBeta=pm.Beta.dist(alpha=tracePoint["alpha"][i],beta=beta-tracePoint["alpha"][i])
        cb[i]=cheatBeta.logp(values[:,i]).eval()
    return np.array(cb).T

def _log_post_trace(trace, model):
    # copy the contents of _log_post_trace function from pymc3/stats.py
    # but replace "var.logp_elemwise(pt)" with "cheat_logp(pt,model)"
    # <...>

def mywaic(trace, model=None, pointwise=False):
    # copy the contents of waic function from pymc3/stats.py
    # <...>

显然，这个骗局不是很“好”，而且我仍然对如何以适当的方式取得类似的结果非常感兴趣。当然，如果可能的话。在