Для нецентрального хи-квадрата в scipy, как определяется вход «nc»?

Question

scipy.stats.ncx2 реализует некоторые функции для нецентрального распределения хи-квадрат. Для этих функций есть один вход 'nc'.

Предположим, есть k независимых случайных чисел из N (mu, 1)

Мой вопрос: следует ли определять nc как k mu^ 2 или \ sqrt (k mu ^ 2).

Я спрашиваю об этом, потому что из Википедии прямо говорится следующее:

"В качестве альтернативы pdf может быть записан как

exp (- (nc + df) / 2) * 1/2 * (x / nc) ** ((df-2) / 4) * I (df-2) / 2

где параметр нецентральности в этой формуле - это квадратный корень из суммы квадратов. "

И в документации scipy.stats.ncx2 pdf имеет точно такую ​​же форму, как указано выше.

Итак, должен ли ввод 'nc' быть суммой квадратов или квадратным корнем из суммы квадрата.

Есть ли какой-либо способ числовой проверки этого?

Answer 1

Значения параметра нецентральности в этих двух представлениях PDF на странице википедии одинаковы. Они не изменили определение λ, которое является суммой квадратов средних нормальных распределений.

Вот скрипт, который генерирует те же кривые, что и график на странице википедии. Цветные линии вычисляются с использованием 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 соответствует нормализованной гистограмме образца.