对于scipy中的非中心chisquare，如何定义输入“nc”？

Question

你知道吗scipy.stats.ncx2实现非中心卡方分布的一些函数。这些函数有一个输入'nc'。你知道吗

假设N（mu，1）中有k个独立的随机数

我的问题是nc应该被定义为kmu^2还是sqrt（kmu^2）。你知道吗

我这样问是因为维基百科明确指出：

“或者，pdf可以写成

实验（-（nc+df）/2）*1/2*（x/nc）**（（df-2）/4）*I（df-2）/2

式中，非中心性参数是平方和的平方根

以及scipy.stats.ncx2，pdf格式与上述格式完全相同。你知道吗

所以输入'nc'应该是平方和，还是平方和的平方根。你知道吗

有什么方法可以从数值上验证这一点吗？你知道吗

Answer 1

wikipedia页面中PDF的这两种表示形式中的noncentrality参数的含义是相同的。他们没有改变λ的定义，λ是正态分布均值的平方和。你知道吗

下面是一个脚本，它生成的曲线与wikipedia页面中的曲线相同。彩色线是用scipy.stats.ncx2.pdf计算的，灰色线是用维基百科页面中给出的无穷级数的前10项计算的。该图验证这些表达式是否只是相同值的不同表达式。你知道吗

import numpy as np
from scipy.stats import ncx2, chi2
import matplotlib.pyplot as plt


def approx_pdf(x, k, lam):
    p = np.zeros_like(x, dtype=np.float64)
    f = 1
    for i in range(10):
        p += np.exp(-lam/2) * (lam/2)**i * chi2.pdf(x, k + 2*i) / f
        f *= (i + 1)
    return p

# df == k on wikipedia
# nc == lambda on wikipedia

x = np.linspace(0, 8, 400)

linestyle = '-'
for df in [2, 4]:
    for nc in [1, 2, 3]:
        plt.plot(x, ncx2.pdf(x, df, nc), linestyle=linestyle,
                 label=f'k = {df}, λ = {nc}')
        plt.plot(x, approx_pdf(x, df, nc), 'k', alpha=0.1, linewidth=6)
    linestyle = ' '

plt.title("Noncentral chi-square distribution\nProbability density function")
plt.xlabel('x')
plt.legend(shadow=True)
plt.grid(alpha=0.3)
plt.tight_layout()
plt.show()

脚本生成的绘图：

---

下面是另一个简短的脚本，来演示非中心性参数，实际上是正态分布平均值的平方和。它生成一个大样本值，每个值是三个正态随机变量的平方和，平均值分别为1、1.5和3。这个样本的分布应该是一个非中心卡方分布，有3个自由度，非中心参数等于平均值的平方和。你知道吗

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


# Means of the normal distributions.
mu = np.array([1, 1.5, 3])

k = len(mu)          # df in scipy.stats.ncx2
lam = (mu**2).sum()  # nc in scipy.stats.ncx2

# The distribution of sample should be a noncentral chi-square
# with len(mu) degrees of freedom and noncentrality sum(mu**2).
sample = (np.random.normal(loc=mu, size=(100000, k))**2).sum(axis=1)

# Plot the normalized histogram of the sample.
plt.hist(sample, bins=60, density=True, alpha=0.4)

# This plot of the PDF should match the histogram.
x = np.linspace(0, sample.max(), 800)
plt.plot(x, ncx2.pdf(x, k, lam))

plt.xlabel('x')
plt.grid(alpha=0.3)
plt.show()

如图中所示，理论PDF与样本的标准化直方图相匹配。你知道吗