标签: probabilistic-programming

我试图用高斯(和更复杂)函数拟合一些数据.我在下面创建了一个小例子.

我的第一个问题是,我做得对吗？

我的第二个问题是,如何在x方向上添加错误,即在观察/数据的x位置？

如何在pyMC中进行这种回归很难找到很好的指南.也许是因为它更容易使用一些最小二乘或类似的方法,但我最终有很多参数,需要看看我们如何约束它们并比较不同的模型,pyMC似乎是一个很好的选择.

import pymc
import numpy as np
import matplotlib.pyplot as plt; plt.ion()

x = np.arange(5,400,10)*1e3

# Parameters for gaussian
amp_true = 0.2
size_true = 1.8
ps_true = 0.1

# Gaussian function
gauss = lambda x,amp,size,ps: amp*np.exp(-1*(np.pi**2/(3600.*180.)*size*x)**2/(4.*np.log(2.)))+ps
f_true = gauss(x=x,amp=amp_true, size=size_true, ps=ps_true )

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

# define the model/function to be fitted.
def model(x, f): 
    amp = pymc.Uniform('amp', …

我使用了基于Clojure的'Anglican',我觉得这对我不好.糟糕的文件和太小的社区找不到帮助.此外,我仍然无法熟悉基于Scheme的语言.所以我想将语言改为基于Python的东西.

也许Pyro或Pymc可能就是这种情况,但我完全不知道这两者.

• 这两个框架有什么区别？
• 他们可以用于同样的问题吗？
• 有没有例子,比较闪耀的是什么？

我有一个n个观察列表,每个观察值是两个Weibull分布式变量的总和:

x[i] = t1[i] + t2[i]
t1[i] ~ Weibull(shape1, scale1)
t2[i] ~ Weibull(shape2, scale2)

我的目标是:

1)估计威布尔分布(shape1,scale1,shape2,scale2)的形状和比例参数,

2)对于每个观测值x [i],估计t1 [i](并且t2 [i]由此得出).

(旁白:每次观察x [i]是癌症诊断的年龄,t1 [i]和t2 [i]是肿瘤发展的两个不同时期.实际模型也涉及突变数据,但在此之前尝试一下,我想确保我可以使用PyMC来解决这个更简单的问题.)

我正在使用PyMC2进行这些估算,看起来运行会收敛,但结果不正确.我不知道我的PyMC模型语法,MCMC设置或两者都有问题.我试着适应这个建议使用电位潜变量模型.首先,我为每个观察定义x [i]和t1 [i]:

for i in xrange(n):
    x[i] = pm.Index('x_%i'%i, x=data, index=i) # data is a list of observations
    t1[i] = pm.Weibull('t1_%i'%i, alpha=shape1, beta=scale1)
    # Ensure that initial guess for t1 is not more than the observed sum:
    if t1[i].value >= x[i].value:
        t1[i].value = 0.95 * x[i].value

然后我为t2 [i] = x [i] …

如何使用pymc参数化概率图形模型？

假设我有两个节点的PGM X和Y.让我们说X->Y是图表.

并且X采用两个值{0,1},并且 Y还采用两个值{0,1}.

我想使用pymc来学习分布的参数,并用它来填充图形模型以进行推理.

我能想到的方式如下:

X_p = pm.Uniform("X_p", 0, 1)
X = pm.Bernoulli("X", X_p, values=X_Vals, observed=True)
Y0_p = pm.Uniform("Y0_p", 0, 1)
Y0 = pm.Bernoulli("Y0", Y0_p, values=Y0Vals, observed=True)
Y1_p = pm.Uniform("Y1_p", 0, 1)
Y1 = pm.Bernoulli("Y1", Y1_p, values=Y1Vals, observed=True)

这里Y0Vals是Y对应于X值= 0 Y1Vals的值,并且是Y对应于X值= 1的值.

计划是从这些中抽取MCMC样本,并使用Y0_p和Y1_p 填充离散贝叶斯网络的概率...所以概率表适用P(X) = (X_p,1-X_p)于P(Y/X):

  Y  0 …

我正在尝试改编Cameron Davidson-Pilon的贝叶斯黑客方法，第1章“介绍我们的第一把锤子：PyMC”中的文本消息示例， 以处理多种观察结果。下面的解决方案似乎正在工作，但是我对pymc并不陌生，我不确定这是否是处理pymc中多个时间序列观测值的好方法。任何建议将不胜感激！

为了重述贝叶斯黑客方法的文本消息示例，观察结果包括74天的文本消息计数，如下图所示。

该书使用一个切换点参数（tau）和两个指数参数（lambda1和lambda2）对该过程进行建模，这两个参数分别控制tau前后的泊松分布消息数。对于此示例，pymc使用以下代码产生大约为tau = 45，lambda1 = 18和lambda2 = 23的解决方案，该代码与本书的代码几乎相同：

import numpy as np
import pymc

observation = np.loadtxt( './txtdata.csv' ) #data available at the book's GitHub site
n_days      = observation.size    #number of days
alpha       = 1./20  #assume a mean of 20 messages per day
lambda1     = pymc.Exponential("lambda1", alpha)
lambda2     = pymc.Exponential("lambda2", alpha)
tau         = pymc.DiscreteUniform("tau", lower=0, upper=n_days)

@pymc.deterministic
def lambda_(tau=tau, lambda1=lambda1, lambda2=lambda2):
    a       = np.zeros(n_days)
    a[:tau] = lambda1
    a[tau:] = lambda2
    return a
observation_model  = pymc.Poisson("observation", …

在这里，我的目的是估算由下式给出的阻尼谐波振荡器的参数（伽马和ω）

dX ^ 2 / dt ^ 2 +γ* dX / dt +（2 * pi * omega）^ 2 * X = 0。（我们可以向系统中添加高斯白噪声。）

 import pymc
 import numpy as np
 import scipy.io as sio
 import matplotlib.pyplot as plt; 
 from scipy.integrate import odeint

 #import data
 xdata = sio.loadmat('T.mat')['T'][0]  #time
 ydata1 = sio.loadmat('V1.mat')['V1'][0]  #  V2=dV1/dt, (X=V1),
 ydata2 = sio.loadmat('V2.mat')['V2'][0]  # dV2/dt=-(2pi*omega)^2*V1-gama*V2

 #time span for solving the equations
 npts= 500   
 dt=0.01
 Tspan=5.0
 time = np.linspace(0,Tspan,npts+1) 

 #initial condition
 V0 = [1.0, 1.0]

# Priors for unknown …

我正在尝试运行以下模型，但在编译过程中失败：

import numpy as np
import pymc3 as pm


def sample_data(G=1, K=2):
    # mean proportion ([0,1]) for each g
    p_g = np.random.beta(2, 2, size=G)

    # concentration around each p_g
    c_g = np.random.lognormal(mean=0.5, sigma=1, size=G)

    # reparameterization for standard Beta(a,b)
    a_g = c_g * p_g / np.sqrt(p_g**2 + (1.-p_g)**2)
    b_g = c_g*(1.-p_g) / np.sqrt(p_g**2 + (1.-p_g)**2)

    # for each p_g, sample K proportions
    p_gk = np.random.beta(a_g[:, np.newaxis], b_g[:, np.newaxis], size=(G, K))

    return p_gk

# Data size
G = 3
K = …

来自 Pytorch-Pyro 的网站：

我们很高兴地宣布发布 NumPyro，这是一个 NumPy 支持的 Pyro，使用 JAX 进行自动微分和 JIT 编译，HMC 和 NUTS 的速度提高了 100 倍以上！

我的问题：

1. NumPyro（超过 Pyro）的性能提升（有时是340 倍或 2 倍）究竟来自哪里？
2. 更重要的是，为什么（而不是，在哪里）我会继续使用 Pyro？

额外的：

1. 与 Tensorflow Probability 相比，我应该如何查看 NumPyro 的性能和特性，以决定在何处使用哪个？