将rJAGS被审查的线性回归模型转换为PyMC3

Question

我正在尝试将我的老板最初用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\u回归）_模型.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中的tilde本身的概念对我来说还是有点奇怪）。你知道吗

如果有人能帮助我，我将不胜感激。据我所知，它甚至可能不可能与PyMC3，或者它很容易，我只是错过了一些东西。不管怎样，我已经把头撞在墙上好几天了，所以我想现在最好还是寻求帮助。你知道吗