问题

I have run a bayesian model (regression) using pymc5.
我已经运行了一个使用pymc5的贝叶斯模型（回归）。

I would like to get a score indicator (AIC|BIC or LOO) to compare different models.
我想获取一个评分指标（AIC|BIC或LOO）来比较不同的模型。

I am a newcomer to pymc and I don't understand how to get this indicator.
我是pymc的新手，不太明白如何获得这个指标。

For reference here the basic code:
以下是基本的代码供参考：

# bayes model
basic_model = pm.Model()

with basic_model:
    
    cost = pm.Normal("cost", mu=0, sigma=2)
    
    a = pm.Normal('a_array', mu=0, sigma=2)
    b = pm.Normal('b_array', mu=0, sigma=2)
    c = pm.HalfNormal('c_array', sigma=2)
    sigma = pm.HalfNormal("sigma", sigma=1)

    # Expected value of outcome
    mu = cost + a * a_array + b * b_array + c * c_array
    # Likelihood (sampling distribution) of observations
    Y_obs = pm.Normal("Y_obs", mu=mu, sigma=sigma, observed=y_array)
    

with basic_model:
    # draw 1000 posterior samples
    idata = pm.sample(1000)

I tried this code:
我尝试了这段代码：

f_logp = basic_model.compile_logp()
# k is the total number of estimated parameters
k = sum(list(map(len, map_est.values())))
aic = 2 * k - 2 * f_logp(map_est)

but it doesn't work (len() of unsized object)
但它不起作用（len()无法用于未定大小的对象）。

Could anyone help me?
有人可以帮助我吗？
All the best
祝一切顺利。

# bayes model
basic_model = pm.Model()

with basic_model:
    
    cost = pm.Normal("cost", mu = 0, sigma=2)
    
    a= pm.Normal('a_array', mu=0, sigma=2)
    b = pm.Normal('b_array', mu=0, sigma=2)
    c= pm.HalfNormal('c_array', sigma=2)
    sigma = pm.HalfNormal("sigma", sigma=1)

    # Expected value of outcome
    mu = cost + a* a_array + b* b_array + c*c_array
    # Likelihood (sampling distribution) of observations
    Y_obs = pm.Normal("Y_obs", mu=mu, sigma=sigma, observed=y_array)
    

with basic_model:
    # draw 1000 posterior samples
    idata = pm.sample(1000)

f_logp = basic_model.compile_logp()
# k is the total number of estimated parameters
k = sum(list(map(len, map_est.values())))
aic  = 2 * k - 2 * f_logp(map_est)

答案1

得分: 2

You can get pymc to compute the log likelihood, and then Arviz for calculating the metrics:

import arviz as az

with basic_model:
    # draw 1000 posterior samples
    idata_01 = pm.sample(1000)

with basic_model:
    pm.compute_log_likelihood(idata_01)

# for loo
basic_model_loo = az.loo(idata_01)

or alternatively passing idata_kwargs={"log_likelihood": True} to sample:

with basic_model:
    # draw 1000 posterior samples
    idata_01 = pm.sample(1000, idata_kwargs={"log_likelihood": True})

# for loo
basic_model_loo = az.loo(idata_01)

See the pymc documentation for more information on model comparison.

import arviz as az

with basic_model:
    # draw 1000 posterior samples
    idata_01 = pm.sample(1000)

with basic_model:
    pm.compute_log_likelihood(idata_01)

# for loo
basic_model_loo = az.loo(idata_01)

with basic_model:
    # draw 1000 posterior samples
    idata_01 = pm.sample(1000, idata_kwargs={"log_likelihood": True})

# for loo
basic_model_loo = az.loo(idata_01)