Неверное поведение со сверткой в ​​PyMC3 / Theano - PullRequest
Новая
       78

Неверное поведение со сверткой в ​​PyMC3 / Theano

0 голосов
/ 07 февраля 2019

В моей первоначальной работе я просто пытался определить подходящее масштабирование для свернутого сигнала и фильтра, который показывает слишком много шума и низкие значения масштабирования.Вывод параметра масштаба работает с линейной регрессией.

Привет!Я новичок в Python, PyMC3 и общем байесовском анализе, поэтому, возможно, я пропустил некоторые основные проблемы, и я был бы очень признателен за помощь.У меня есть электронная система, которую я пытаюсь охарактеризовать с помощью одного фильтра.Входной сигнал свернут с фильтром, чтобы получить измеренный выход.У меня есть отдельный, ручной анализ, который дает результат, но он менее чем удовлетворительный, и я хочу попробовать байесовские подходы.

Я хотел убедиться, что я понимаю различные шаги, прежде чем погружаться дальше, поэтому ясделал фиктивный тест-кейс.Я создал произвольную функцию фильтра, которая свернута с входным сигналом.Чтобы ввести некоторые переменные логического вывода, я умножаю результат на параметр масштабирования, а затем моя функция правдоподобия представляет собой нормальное распределение для свертанного сигнала со стандартным отклонением, которое сравнивается с наблюдением.Поскольку я создал данные, я знаю, какие ответы должны быть получены.Я могу сделать тот же процесс с линейной регрессией, и он правильно получает правильное масштабирование и стандартное отклонение для шума. 

import numpy as np
import matplotlib.pyplot as plt
import pymc3 as pm
import pandas as pd
import theano.tensor as tt
import theano
import theano.tensor.signal.conv

#theano.config.compute_test_value = 'off'

#create source function
dtbase = 0.5e-9
tbase = np.arange(0, 1e-7, dtbase)

set1 = ((tbase > 5e-9) & (tbase < 45e-9))
sig1 = np.zeros((len(tbase),1))
sig1[set1,0] = 1

set2 = ((tbase > 5e-9) & (tbase < 8e-9))
sig2 = np.zeros((len(tbase),1))
sig2[set2,0] = 1

#define the time base of the transfer function
t2 = np.arange(0, 50e-9, dtbase)
hfunc = np.zeros((len(t2),1))

mu1 = 2e-9
sd1 = 0.5e-9
g1 = np.exp(-0.5*(t2 - mu1)**2/sd1**2)

RC1 = 2e-9
A1 = 1
mu2 = 2e-9
RC2 = 20e-9
A2 = 0.02

set3 = (t2 > mu2)
hfunc[set3,0]= A1*np.exp(-1*t2[set3]/RC1) + A2*np.exp(-1*t2[set3]/RC2)
hfunc2 = np.convolve(hfunc[:,0], g1)
hfunc2 = hfunc2[0:len(t2)]/np.trapz(hfunc2)

hproto = np.zeros((len(hfunc2),1))
hproto[:,0] = hfunc2
hfunc3 = hfunc2/(max(hfunc2))

#create the observed data arrays
obs1 = np.convolve(sig1[:,0], hfunc2)
obs1 = obs1[0:len(tbase)]
noiseAmp = np.random.normal(0, 0.005*max(obs1), len(tbase))
obs1 += noiseAmp

obs2 = np.convolve(sig2[:,0], hfunc2)
obs2 = obs2[0:len(tbase)]
noiseAmp2 = np.random.normal(0, 0.005*max(obs1), len(tbase))
obs2 += noiseAmp2



#define the theano convolution function
conv2d = tt.signal.conv.conv2d

x = tt.dmatrix('x')
x.tag.test_value = np.random.rand(5,1)

y = tt.dmatrix('y')
y.tag.test_value = np.random.rand(5,1)
veclen = x.shape[0]

conv1d_expr = conv2d(x, y, \
                 image_shape=(1, veclen), \
                 border_mode='full')

conv1d = theano.function([x, y], \
                     outputs=conv1d_expr)

#append the two signals to process at the same time
xapp = np.append(sig1, sig2)
obsapp = np.append(obs1, obs2)

basic_model = pm.Model()

with basic_model:

    scale = pm.Normal('scale', mu=0, sd = 5)

    Yobs_a = conv1d(sig1, hproto)
    Yobs_b = conv1d(sig2, hproto)

    #Yobs = (Yobs_a[0:len(sig1)])*(1-scale)
    Ydelta_a = Yobs_a[0:len(sig1)]*(1-scale)
    Ydelta_b = Yobs_b[0:len(sig2)]*(1-scale)

    #Ydelta = th_append(Ydelta_a, Ydelta_b)    
    Ydelta = tt.concatenate([Ydelta_a, Ydelta_b], axis=0)

    sigma = pm.HalfCauchy('sigma', 1)

    obs = pm.Normal('obs', mu=Ydelta, sd=sigma, observed=obsapp)

    trace = pm.sample(800, tune=800, cores=1, jobs=1, chains=2)

Теперь, если я использую: 

ppc = pm.sample_posterior_predictive(trace, samples=100, model=basic_model)

for i in range(100):
    plt.plot((ppc['obs'][i,:]), '-b', alpha=0.05)
    plt.plot(obsapp, '-r', linewidth=1)

Я получаюочень шумный набор постеров, потому что мой решатель возвращается со значением: 

       mean   sd  mc_error  hpd_2.5  hpd_97.5    n_eff  Rhat
scale  0.72  0.0       0.0     0.72      0.73  1549.32   1.0
sigma  0.38  0.0       0.0     0.38      0.38  1641.81   1.0

, тогда как масштаб должен быть 0, а сигма должна быть 0,005!

Если я запускаю тот же алгоритм, используямасштабированная линейная регрессия как «наблюдаемый», этот метод возвращается с очень точными значениями.

Я явно что-то упускаю, но не знаю, где искать.

...