使用Crystal Ball函数时的异常绘图问题：曲线存在不连续性

你的代码出现不连续的核心原因是高斯部分和幂律部分的缩放因子不匹配：Crystal Ball函数的高斯部分应基于标准化后的变量（( z=(x-\mu)/\sigma )）使用标准正态分布的概率密度，而你直接用了norm.pdf(x, loc=mean, scale=sigma)——这个函数会自动除以(\sigma)，导致高斯部分比幂律部分多了一个(1/\sigma)的缩放因子，最终在衔接点（( z=-\alpha )）处数值不连续。

修正方案 #

有两种等价的修正方式，任选其一即可：

方式1：调整高斯部分为标准正态pdf（基于z值） #

将高斯部分的计算改为基于标准化变量z，直接调用标准正态分布的pdf（即norm.pdf(z)，默认loc=0, scale=1）：

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

def crystal_ball(x, alpha, n, mean, sigma, N):
    alpha = abs(alpha)
    z = (x - mean) / sigma
    A = (n / alpha) ** n * np.exp(-0.5 * alpha ** 2)
    B = n / alpha - alpha

    mask = z > -alpha

    # 高斯部分：用标准正态pdf（基于z）
    gaussian_part = N * norm.pdf(z[mask])
    # 幂律部分：保持原有计算
    powerlaw_part = N * A / ((B - z[~mask]) ** n)

    result = np.zeros_like(x)
    result[mask] = gaussian_part
    result[~mask] = powerlaw_part
    return result, mask

# 示例调用
x_values = np.linspace(-10, 4, 1000)
alpha_val = 1
n_val = 10
mean_val = 0
sigma_val = 1
N_val = 1

y_values, mask = crystal_ball(x_values, alpha_val, n_val, mean_val, sigma_val, N_val)

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='powerlaw')
plt.vlines(-alpha_val, 0, 1, colors='orange', lw=1)
plt.xlabel('Invariant Mass')
plt.ylabel('Probability Density')
plt.legend()
plt.show()

方式2：给幂律部分添加1/σ的缩放因子 #

如果想保留norm.pdf(x, loc=mean, scale=sigma)的写法，只需给幂律部分除以(\sigma)，匹配高斯部分的缩放：

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

def crystal_ball(x, alpha, n, mean, sigma, N):
    alpha = abs(alpha)
    z = (x - mean) / sigma
    A = (n / alpha) ** n * np.exp(-0.5 * alpha ** 2)
    B = n / alpha - alpha

    mask = z > -alpha

    # 高斯部分：保留原有写法
    gaussian_part = N * norm.pdf(x[mask], loc=mean, scale=sigma)
    # 幂律部分：添加1/sigma的缩放
    powerlaw_part = N * A / ((B - z[~mask]) ** n) / sigma

    result = np.zeros_like(x)
    result[mask] = gaussian_part
    result[~mask] = powerlaw_part
    return result, mask

# 示例调用同方式1

验证效果 #

用你提供的示例参数运行修正后的代码，会得到完全连续的Crystal Ball曲线，高斯部分和幂律部分在( z=-\alpha )（即x=-1）处完美衔接。

内容的提问来源于stack exchange，提问作者Magic Dealer