如何编写一个基本的 Pymc3 模型？教程中的一个小替换是有问题的

问题描述

我正在尝试在给定确定性方程和模拟噪声数据的情况下推断 2 个参数（beta 和 gamma）。出于某种原因，我使用的方程似乎有问题，因为我只是复制了基本的 pymc3 教程并使用了我自己的确定性方程。这是我正在使用的模型：

# True parameter values
beta,gamma = 0.21,0.07

# Size of dataset
days = 50

# Predictor variable
time = np.arange(0,days,1)

# Simulate outcome variable
data = []
for t in time:
    data.append((beta/((beta-gamma))*(np.exp(t*(beta-gamma))-1)+1) + np.random.normal(0,1))

basic_model = pm.Model()

def smodel(beta,gamma):
    s = beta/((beta-gamma))*(tt.exp(time*(beta-gamma))-1)+1
    return s

with basic_model:

    # Priors for unkNown model parameters
    beta = pm.normal("beta",mu=0,sigma=10)
    gamma = pm.normal("gamma",sigma=10)

    # Expected value of outcome
    #smodel_pm = pm.Deterministic('smodel',smodel(inputParam))
    y_obs = pm.normal('obs',mu=smodel(beta,gamma),sigma=1,observed=data)
    
    # Draw the specified number of samples
    trace = pm.sample(step=pm.Metropolis())

但是，当我运行跟踪摘要时，所有内容都为 0。有人知道是什么问题吗？

解决方法

这里的参数化相当差（相关变量，域中的对称解），加上 Metropolis-Hastings 只需要运行很长时间，而默认设置假设为 NUTS。

以下是建议的替代参数化，以及对于此采样策略更合理的调整和绘制计数：

basic_model = pm.Model()

def smodel(a,b):
    s = a*(tt.exp(b*time)-1)+1
    return s

with basic_model:

    # priors for pre-transformed model parameters
    a = pm.Normal("a",mu=0,sigma=10)
    b = pm.HalfNormal("b",sigma=10)
    
    # (transformed) parameters of interest
    beta  = pm.Deterministic("beta",a*b)
    gamma = pm.Deterministic("gamma",(a-1)*b)
    
    # expected value of outcome
    y_obs = pm.Normal('obs',mu=smodel(a,b),sigma=1,observed=data)
    
    # Draw the specified number of samples
    trace = pm.sample(step=pm.Metropolis(),tune=100000,draws=50000)

人们可能会进一步提高 draws，因为有效样本量 (ESS) 太小了，但感兴趣的参数的数字应该是：

查看配对图，可以看到 a,b 的相关性仍然很高，这说明样本高度自相关。

痕迹和密度（每条链）对我来说看起来不错：