使用水晶球功能时出现问题 - 意外的绘图行为

Question

我正在尝试使用水晶球函数来拟合热量计数据中的一些能谱。由于这是我第一次使用它，我之前尝试过一些测试。我在维基百科（https://en.wikipedia.org/wiki/Crystal_Ball_function）和其他文献中查阅了水晶球函数的定义，然后通过以下Python代码实现了它：

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm

def crystal_ball(x, alpha, n, mean, sigma, N):

    A = (n / abs(alpha)) ** n * np.exp(-0.5 * alpha ** 2)
    B =  n / abs(alpha) - abs(alpha)

    mask = ((x - mean) / sigma) > -alpha

    gaussian_part = N * norm.pdf(x[mask], loc = mean, scale = sigma)
    powerlaw_part = N * A / ((B - (x[~mask] - mean) / sigma) ** n)

    result = np.zeros_like(x)

    result[mask]  = gaussian_part
    result[~mask] = powerlaw_part        

    return result, mask

# Example usage:
x_values = np.linspace(-10, 4, 1000)
alpha_val = 1
n_val = 10
mean_val = 0
sigma_val = 1
N_val = 1

y_values = crystal_ball(x_values, alpha_val, n_val, mean_val, sigma_val, N_val)[0]

mask = crystal_ball(x_values, alpha_val, n_val, mean_val, sigma_val, N_val)[1]

plt.plot(x_values, y_values, c = 'b', marker = '*', ms = 1, ls = '',
         label = f'Crystal Ball (α={alpha_val}, n={n_val}, μ={mean_val}, σ={sigma_val})')

plt.plot(x_values[mask], y_values[mask], c = 'g', marker = '*', ms = 1, ls = '',
         label = 'gauss')

plt.plot(x_values[~mask], y_values[~mask], c = 'r', marker = '*', ms = 1, ls = '',
         label = 'exp')

plt.vlines(-alpha_val, 0, 1, colors = 'orange', lw = 1)

plt.xlabel('Invariant Mass')
plt.ylabel('Probability Density')
plt.legend()
plt.show()

但是在绘图时我遇到了意想不到的不连续性：

就好像函数的幂律部分乘以一个我既无法找到也无法解释的比例因子。有人可以解释这种行为的原因，或者至少告诉我代码中的错误吗？

在实践中，我希望获得连续的曲线。

我也尝试过使用

scipy.stats.crystalball

模块，但结果发现它更加复杂并且没有用。

Answer 1

你的问题是你创建了不一致的值。根据维基百科页面，您只能控制

alpha

、

n

、

xbar

和

sigma

的值，但当该值实际上应由其他值确定时，您也可以设置

N

。如果您删除

N

定义并生成每个维基百科页面的所有值，那么您将获得正确的结果。

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import erf

plt.close("all")

def crystal_ball(x, alpha, n, xbar, sigma):
    n_over_alpha = n/abs(alpha)
    exp = np.exp(-0.5*alpha ** 2)
    A = (n_over_alpha)**n*exp
    B =  n_over_alpha - abs(alpha)
    C = n_over_alpha/(n-1)*exp
    D = np.sqrt(0.5*np.pi)*(1 + erf(abs(alpha)/np.sqrt(2)))
    N = 1/(sigma*(C + D))

    mask = (x - xbar)/sigma > -alpha

    gaussian = N*np.exp(-0.5*((x[mask]-xbar)/sigma)**2)
    powerlaw = N*A*(B - (x[~mask]-xbar)/sigma)**-n

    result = np.zeros_like(x)
    result[mask] = gaussian
    result[~mask] = powerlaw
    # or
    # result = np.concatenate((powerlaw, gaussian))

    return result, mask

x = np.linspace(-10, 4, 1000)
alpha = 1
n = 10
xbar = 0
sigma = 1

y, mask = crystal_ball(x, alpha, n, xbar, sigma)

plt.plot(x[mask], y[mask], label="gauss")
plt.plot(x[~mask], y[~mask], label="exp")
plt.legend()

现在，您似乎夸大了使用 scipy 的

crystallball

函数的难度。如果您查看文档，您会发现他们所说的

beta

就是您所说的

alpha

，他们所说的

m

就是您所说的

n

。文档都说，

上面的概率密度是以“标准化”形式定义的。要移动和/或缩放分布，请使用

loc

和

scale

参数。具体来说，

crystalball.pdf(x, beta, m, loc, scale)

与

crystalball.pdf(y, beta, m) / scale

与

y = (x - loc) / scale

完全相同。 换句话说，他们已经完全实现了您正在寻找的内容，只是名称不同（

loc

是

xbar

，

scale

是

sigma

）。

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import erf
from scipy import stats

plt.close("all")

def crystal_ball(x, alpha, n, xbar, sigma):
    n_over_alpha = n/abs(alpha)
    exp = np.exp(-0.5*alpha ** 2)
    A = (n_over_alpha)**n*exp
    B =  n_over_alpha - abs(alpha)
    C = n_over_alpha/(n-1)*exp
    D = np.sqrt(0.5*np.pi)*(1 + erf(abs(alpha)/np.sqrt(2)))
    N = 1/(sigma*(C + D))

    mask = (x - xbar)/sigma > -alpha

    gaussian = N*np.exp(-0.5*((x[mask]-xbar)/sigma)**2)
    powerlaw = N*A*(B - (x[~mask]-xbar)/sigma)**-n

    result = np.zeros_like(x)
    result[mask] = gaussian
    result[~mask] = powerlaw
    # or
    # result = np.concatenate((powerlaw, gaussian))

    return result, mask

x = np.linspace(-10, 4, 100)
alpha = 1
n = 10
xbar = 0
sigma = 1

y, mask = crystal_ball(x, alpha, n, xbar, sigma)
y_scipy = stats.crystalball.pdf(x, beta=alpha, m=n, loc=xbar, scale=sigma)

plt.plot(x[mask], y[mask], "x", label="gauss")
plt.plot(x[~mask], y[~mask], "x", label="exp")
plt.plot(x, y_scipy, label="scipy")
plt.legend()

如您所见，结果是相同的。