pymc3中非对称高斯参数的拟合

Question

我正在尝试将非对称高斯分布拟合到我的数据中。我的数据只是一个名为wave（x）的numpy数组和一个名为spec（y）的numpy数组，看起来像一个不对称的高斯分布

This is the image with the data with an asymmetric gaussian fitted with curve_fit (this has a continuum too, but this is not important right now.

这就是功能：

def agauss(amp, cen, b_sigma, r_sigma, x):
y = np.zeros(len(x))
for i in range(len(x)):
    if x[i] < cen:
        y[i] = amp*np.exp(-((x[i] - cen)**2)/(2*b_sigma**2))
    else:
        y[i] = amp*np.exp(-((x[i] - cen)**2)/(2*r_sigma**2))
return y

我使用此代码来拟合参数：

with pm.Model() as asym:
    cen = pm.Uniform('cen', lower=5173, upper=5179)
    bsigma = pm.HalfCauchy('bsigma', beta=3)
    rsigma = pm.HalfCauchy('rsigma', beta=3)
    amp = pm.Uniform('amp', lower=1e-19, upper=1e-16)

    err = pm.HalfCauchy('err', beta=0.0000001)

    ag_pred = pm.Normal('ag_pred', mu=agauss(amp, cen, bsigma, rsigma, wave), sigma=err, observed=spec) 

    agdata = pm.sample(3000, cores=2)

但是我在theano.tensor模块中得到错误“变量不支持布尔运算”。我应该如何定义函数以适应参数？有更好的办法吗？谢谢

    144     err = pm.HalfCauchy('err', beta=0.0000001)
    145 
--> 146     ag_pred = pm.Normal('ag_pred', mu=agauss(amp, cen, bsigma, rsigma, wave), sigma=err, observed=spec)
    147 
    148     agdata = pm.sample(3000, cores=2)

~/Documents/OIII_emitters/m2fs_reduction/test/assets/scripts/analysis.py in agauss(amp, cen, b_sigma, r_sigma, x)
    118     y = np.zeros(len(x))
    119     for i in range(len(x)):
--> 120         if x[i] < cen:
    121             y[i] = amp*np.exp(-((x[i] - cen)**2)/(2*b_sigma**2))
    122         else:

~/anaconda3/envs/data_science/lib/python3.7/site-packages/theano/tensor/var.py in __bool__(self)
     92         else:
     93             raise TypeError(
---> 94                 "Variables do not support boolean operations."
     95             )
     96 

TypeError: Variables do not support boolean operations.

Answer 1

下面是一种遵循pymc3的switchpoint技巧的方法。这基本上是thisSO问题的非线性扩展。下面我展示了代码，以及它如何处理我创建的一些模拟数据。我不希望模拟数据与实际数据太相似，但它应该为您提供一个起点。这种方法确实存在一些建模限制（即拟合两个高斯），但同样，它们可能会随着更真实的输入数据而消失

首先，我创建一些模拟输入数据。请注意，相互连接的两个高斯将具有不同的法线，并且在某种程度上也将以不同的方式接近光谱的连续性。我没有试图在模拟数据中包含后一个细节，但我保持了标准化的一致性。另一个复杂因素是两个高斯的振幅不同。这对于实际数据可能不会有问题，但在这里我只是添加一个常量来匹配它们。我们希望从这些模拟数据恢复参数cen_mock、sigma_mock_1和sigma_mock_2。最后请注意，我保持了与您问题中相同的频率限制

import numpy as np
import pymc3 as pm
import theano.tensor as tt
import matplotlib.pyplot as plt


def gaussian_pdf(x,sigma,cen):
    g1 = (1/(sigma*np.sqrt(2*np.pi)))*np.exp(
         -0.5*np.square(x-cen)/np.square(sigma))
    return g1

# create mock data
wavelength_min = 5173
wavelength_max = 5179
cen_mock = 5175
x1 = np.arange(5173,cen_mock,0.01)
x2 = np.arange(cen_mock,5179,0.01)
x = np.arange(wavelength_min,wavelength_max,0.01)

sigma_mock_1 = 1.4
g1 = gaussian_pdf(x[x<=cen_mock],sigma_mock_1,cen_mock)
sigma_mock_2 = 1.9
g2 = gaussian_pdf(x[x>cen_mock],sigma_mock_2,cen_mock)

一旦我们有了一些模拟数据，我们就在pymc3和sample中构建模型：

 # construct model
with pm.Model() as asym:
    switchpoint = pm.DiscreteUniform("switchpoint",lower=x.min(),upper=x.max())
    sigma = pm.HalfCauchy('sigma', beta=10, testval=1.)

    sigma_1 = pm.HalfNormal('sigma_1', sd=10)
    sigma_2 = pm.HalfNormal('sigma_2', sd=10)

    sd = pm.math.switch(switchpoint > x, sigma_1, sigma_2)

    likelihood = pm.Normal(
        'y', mu=(1/(sd*np.sqrt(2*np.pi)))*tt.exp(
        -0.5*tt.square(x-switchpoint)/tt.square(sd)), 
        sd=sigma, observed=y_mock)

with asym:
    step1 = pm.NUTS([sigma_1,sigma_2])
    step2 = pm.Metropolis([switchpoint])
    trace = pm.sample(4000, step=[step1,step2],cores=4,tune=4000,chains=4)

模型的参数是swithcpoint，它表示高斯变化的位置和标准偏差。因此，采样器将根据x的值估计“右”或“左”标准偏差

我们现在检查结果：

df = pm.summary(trace)
x1_m = x[x<=df.loc['switchpoint']['mean']]
x2_m = x[x>df.loc['switchpoint']['mean']]
g1_m = gaussian_pdf(x1_m,df.loc['sigma_1']['mean'],df.loc['switchpoint']['mean'])
g2_m = gaussian_pdf(x2_m,df.loc['sigma_2']['mean'],df.loc['switchpoint']['mean'])

amp_m = g1_m.max() - g2_m.max()
plt.plot(x[x<=cen_mock],g1,label='left gaussian mock')
plt.plot(x[x>cen_mock],g2,label='right gaussian mock')
plt.plot(x1_m,g1_m,label='left gaussian model')
plt.plot(x2_m,g2_m+amp_m,label='right gaussian model')
plt.text(5177,0.23,'sd_mock1='+str(sigma_mock_1))
plt.text(5177,0.2,'sd_mock2='+str(sigma_mock_2))
plt.legend(frameon=False)

当然有更好的方法来检查抽样结果（实际上是贝叶斯的），但这是一种快速而肮脏的方法来了解正在发生的事情。请注意，在建模之后，我仍然需要在其中一个高斯数中添加一个常数，以自然地加入它们。对于三对不同的标准偏差，我得到以下曲线图：

现在进行一些讨论。为了评估这种方法有多好，我们需要更好地了解真实数据。正如您所看到的，如果标准偏差显著不同，则模型开始失败。然而，对于相似的标准偏差和标准偏差相等的边缘情况，模型在某种程度上是可以接受的。此外，另一个重要输入是两个高斯数的频率数据点数量是否相等。这将决定每个高斯函数如何接近连续统。它还将确定是否需要在第二个高斯函数中任意添加一个常数，以便在拟合真实数据时以更自然的方式将它们连接起来

总之，这是一种非常符合您期望的参数化的方法，但是在使用模拟数据进行评估时需要做一些额外的工作。从您的问题来看，似乎需要switchpoint方法，但只有当应用于真实数据（或更真实的模拟数据）时，才能确定它是否足够