Использование PyMC3 для растянутой экспоненты: плохая начальная энергия - PullRequest
0 голосов
/ 27 июня 2019

Я пытаюсь изменить самое простое начало работы - пример 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)

Вместо растянутой экспоненты я также попробовал степенные законы и функции синуса. Мне кажется, что проблема возникает, как только моя модель не является инъективной. Может ли это быть проблемой (как видно, я новичок в этой области)? Могу ли я ограничить выборку только положительными значениями х? Есть ли какие-нибудь хитрости к этому?

1 Ответ

0 голосов
/ 28 июня 2019

Таким образом, проблема в том, что

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 был положительным, но мне нужно больше подробностей, чтобы помочь вам описать модель.

...