使用 PyMC3 拟合拉伸指数：不良初始能量

Question

我正在尝试改变最简单的入门 - pymc3 (https://docs.pymc.io/notebooks/getting_started.html) 的示例，这是线性回归拟合拉伸指数的激励示例。

我试过的模型的最简单版本是 y = exp(-x**beta)

import numpy as np
import matplotlib.pyplot as plt
plt.style.use('seaborn-darkgrid')

# Initialize random number generator
np.random.seed(1234)

# True parameter values
sigma = .1
beta = 1

# Size of dataset
size = 1000

# Predictor variable
X1 = np.random.randn(size)

# Simulate outcome variable
Y = np.exp(-X1**beta) + np.random.randn(size)*sigma

# specify the model
import pymc3 as pm
import theano.tensor as tt
print('Running on PyMC3 v{}'.format(pm.__version__))

basic_model = pm.Model()

with basic_model:

    # Priors for unknown model parameters
    beta = pm.HalfNormal('beta', sigma=1)
    sigma = pm.HalfNormal('sigma', sigma=1)

    # Expected value of outcome
    mu = pm.math.exp(-X1**beta)

    # Likelihood (sampling distribution) of observations
    Y_obs = pm.Normal('Y_obs', mu=mu, sigma=sigma, observed=Y)

with basic_model:
    # draw 500 posterior samples
    trace = pm.sample(500)

产生输出

Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [sigma, beta]
Sampling 4 chains:   0%|          | 0/4000 [00:00<?, ?draws/s]/opt/conda/lib/python3.7/site-packages/numpy/core/fromnumeric.py:2920: RuntimeWarning: Mean of empty slice.
  out=out, **kwargs)
/opt/conda/lib/python3.7/site-packages/numpy/core/fromnumeric.py:2920: RuntimeWarning: Mean of empty slice.
  out=out, **kwargs)

Bad initial energy, check any log probabilities that are inf or -inf, nan or very small:
Y_obs   NaN
---------------------------------------------------------------------------
RemoteTraceback                           Traceback (most recent call last)
RemoteTraceback: 
"""
Traceback (most recent call last):
  File "/opt/conda/lib/python3.7/site-packages/pymc3/parallel_sampling.py", line 160, in _start_loop
    point, stats = self._compute_point()
  File "/opt/conda/lib/python3.7/site-packages/pymc3/parallel_sampling.py", line 191, in _compute_point
    point, stats = self._step_method.step(self._point)
  File "/opt/conda/lib/python3.7/site-packages/pymc3/step_methods/arraystep.py", line 247, in step
    apoint, stats = self.astep(array)
  File "/opt/conda/lib/python3.7/site-packages/pymc3/step_methods/hmc/base_hmc.py", line 144, in astep
    raise SamplingError("Bad initial energy")
pymc3.exceptions.SamplingError: Bad initial energy
"""

The above exception was the direct cause of the following exception:

SamplingError                             Traceback (most recent call last)
SamplingError: Bad initial energy

The above exception was the direct cause of the following exception:

ParallelSamplingError                     Traceback (most recent call last)
<ipython-input-310-782c941fbda8> in <module>
      1 with basic_model:
      2     # draw 500 posterior samples
----> 3     trace = pm.sample(500)

/opt/conda/lib/python3.7/site-packages/pymc3/sampling.py in sample(draws, step, init, n_init, start, trace, chain_idx, chains, cores, tune, progressbar, model, random_seed, discard_tuned_samples, compute_convergence_checks, **kwargs)
    435             _print_step_hierarchy(step)
    436             try:
--> 437                 trace = _mp_sample(**sample_args)
    438             except pickle.PickleError:
    439                 _log.warning("Could not pickle model, sampling singlethreaded.")

/opt/conda/lib/python3.7/site-packages/pymc3/sampling.py in _mp_sample(draws, tune, step, chains, cores, chain, random_seed, start, progressbar, trace, model, **kwargs)
    967         try:
    968             with sampler:
--> 969                 for draw in sampler:
    970                     trace = traces[draw.chain - chain]
    971                     if (trace.supports_sampler_stats

/opt/conda/lib/python3.7/site-packages/pymc3/parallel_sampling.py in __iter__(self)
    391 
    392         while self._active:
--> 393             draw = ProcessAdapter.recv_draw(self._active)
    394             proc, is_last, draw, tuning, stats, warns = draw
    395             if self._progress is not None:

/opt/conda/lib/python3.7/site-packages/pymc3/parallel_sampling.py in recv_draw(processes, timeout)
    295             else:
    296                 error = RuntimeError("Chain %s failed." % proc.chain)
--> 297             raise error from old_error
    298         elif msg[0] == "writing_done":
    299             proc._readable = True

ParallelSamplingError: Bad initial energy

INFO (theano.gof.compilelock): Waiting for existing lock by process '30255' (I am process '30252')
INFO (theano.gof.compilelock): To manually release the lock, delete /home/jovyan/.theano/compiledir_Linux-4.4--generic-x86_64-with-debian-buster-sid-x86_64-3.7.3-64/lock_dir
/opt/conda/lib/python3.7/site-packages/numpy/core/fromnumeric.py:2920: RuntimeWarning: Mean of empty slice.
  out=out, **kwargs)
/opt/conda/lib/python3.7/site-packages/numpy/core/fromnumeric.py:2920: RuntimeWarning: Mean of empty slice.
  out=out, **kwargs)

除了拉伸指数，我还尝试了幂律和正弦函数。在我看来，只要我的模型不是单射的，问题就会出现。这会成为一个问题吗（很明显，我是这个领域的新手）？我可以将采样限制为仅正 x 值吗？这有什么技巧吗？

Answer 1

所以这里的问题是

X1**beta

仅在X1 >= 0或beta为整数时定义。当你将其输入到你的观察中时，对于大多数地方，beta 将是一个浮点数，而且很多

mu = pm.math.exp(-X1**beta)

将是 nan。

我用

发现了这个

>>> basic_model.check_test_point()

beta_log__    -0.77
sigma_log__   -0.77
Y_obs           NaN
Name: Log-probability of test_point, dtype: float64

我不确定您要指定的型号！有多种方法可以要求 beta 为整数，也有多种方法可以要求 X1 为正数，但我需要更多详细信息来帮助您描述模型。

使用 PyMC3 拟合拉伸指数：不良初始能量

Using PyMC3 to fit a stretched exponential: bad initial energy

python

regression

pymc3