1

我目前正在尝试将我老板最初用 R/rJAGS 编写的程序翻译成 Python/PyMC3，部分是因为他想看看这是否是 python 可以做的事情，部分是因为我想学习如何做这种事情，知道这似乎是件好事。我已经在 PyMC3 中得到了一个线性拟合模型，但我在尝试复制审查位时遇到了困难。

R 程序读入一个表格，每行具有三个 y 值，用于三个特定的 x 值，这些 x 值在整个数据集中是恒定的。每个 y 值也有一些与之相关的误差。如果是这样，那么我有一个 PyMC3 模型可以做到这一点；这是我为它设置的玩具模型：

import numpy as np
import pymc3 as pmc

# set random seed for reproducibility
np.random.seed(12345)

x = np.linspace(0,10,3)

# Make some model data
# Parameters for linear fit
slope_true = -0.2
inter_true = 0.1

#Linear function
linear = lambda x,slope,inter: slope*x+inter
f_true = linear(x=x,slope=slope_true, inter=inter_true )

# add noise to the data points
f = f_true + np.random.normal(size=len(x)) * 0.05 
f_error = np.ones_like(f_true)*f.max()*np.random.uniform(0,1,size=len(x))

with pmc.Model() as model3:
    slope = pmc.Normal('slope', mu=0, tau=0.4, testval= 0.15)
    inter = pmc.Normal('inter', mu=0, tau=40, testval=0.15)

    linear = pmc.Deterministic('linear', slope*x+inter)

    y = pmc.Normal('y', mu=linear, tau=1.0/f_error**2, observed=f)

    start = pmc.find_MAP()
    step = pmc.NUTS()
    trace = pmc.sample(1000,start=start)

# extract results
slope_fit = np.median(trace.slope)
slope_up = slope_fit - np.percentile(trace.slope, 15.9)
slope_dn = np.percentile(trace.slope, 84.1) - slope_fit

上面的模型是从我在网上找到的例子中拼凑起来的，它在一条线上生成点，添加一些噪音和一些“错误”，然后对有错误的噪音点进行拟合。之后，它会获取斜率的中值以及与之相关的一些错误。

但现在我需要能够解释这些有时会弹出的审查点。在这种情况下，某些 y 值可能未被检测到，因此该点的值被视为审查限制，然后将该点设置为 NaN，但仍与该点相关联的错误。处理此问题的 R 代码模型（另存为 lin_regress_model.bug）如下所示：

model {
    for (i in 1:N) {
        isCensored[i] ~ dinterval(rv[i], censorLimitVec[i])
        rv[i] ~ dnorm(y[i],rve[i])
        y[i] <- a*x[i] + b
    }
    a ~ dnorm(0, 1e-6)
    b ~ dnorm(0, 1e-6)
    tau ~ dgamma(0.001, 0.001)
    sigma <- 1/sqrt(tau)
}

以下是它可能会获取的数据示例：

N = 3 # always 3, because 3 points
isCensored = c(False, False, True)
censorLimitVec = c(-6.65, -6.65, -6.65) # was value of 3rd point before NA
rv = c(-3.4, -4.7, NA) # y-values
rve = c(7e3, 7e2, 6.66) # these are Tau I think, like 1/sigma^2
x = c(0.15, 0.68, 0.94) # x-values

所以所有这些都被传递到了 jags 模型中，它能够适应这些被审查的数据，但我终其一生都无法弄清楚如何将该位转换为 PyMC3 语言。听起来这里的 dinterval 函数可能类似于 PyMC3 中的 Uniform，但我真的不知道该怎么做，因为我不能直接翻译公式行（R 中波浪号本身的概念仍然是对我来说有点奇怪）。

如果有人可以帮助我，将不胜感激。据我所知，PyMC3 甚至可能无法实现，或者它很容易，我只是错过了一些东西。无论如何，我已经把头撞在墙上好几天了，所以我认为此时最好寻求帮助。

4

0 回答 0