0

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

这是带有拟合曲线拟合的不对称高斯数据的图像（这也有一个连续体，但这现在并不重要。

这是功能：

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.
4

1 回答 1

0

这是一种遵循switchpointpymc3 技巧的方法。它基本上是这个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方法是必需的，但只有在应用于真实数据（或更现实的模拟数据）时，才能确定是否足够。

于 2021-07-25T08:24:44.360 回答