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

如何解决在pymc3中拟合非对称高斯参数

我正在尝试将非对称高斯拟合到我的数据中。我的数据只是一个叫做 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,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',observed=spec)
    147 
    148     agdata = pm.sample(3000,cores=2)

~/Documents/OIII_emitters/m2fs_reduction/test/assets/scripts/analysis.py in agauss(amp,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.

解决方法

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

首先我创建了一些模拟输入数据。请注意，相互连接的两个高斯将具有不同的归一化，并且也会以不同的方式接近您的频谱连续体。我没有尝试在模拟数据中包含后者的微妙之处，但我保持规范化一致。另一个复杂因素是两个高斯分布的幅度不同。这可能不会成为真实数据的问题，但在这里我只是添加一个常量来匹配它们。我们希望从这些模拟数据中恢复参数 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 中构建模型并进行采样：

 # 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 方法是必需的，但只有在应用于真实数据（或更真实的模拟数据）时，人们才能确定它是否足够。