В моей первоначальной работе я просто пытался определить подходящее масштабирование для свернутого сигнала и фильтра, который показывает слишком много шума и низкие значения масштабирования.Вывод параметра масштаба работает с линейной регрессией.
Привет!Я новичок в 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!
Если я запускаю тот же алгоритм, используямасштабированная линейная регрессия как «наблюдаемый», этот метод возвращается с очень точными значениями.
Я явно что-то упускаю, но не знаю, где искать.