#python #probability #gaussian
#python #вероятность #гауссовский
Вопрос:
Кажется, я не могу понять, почему теоретические и смоделированные результаты настолько различаются для распределения вероятности нормальной случайной величины в квадрате. (например, мощность сигнала напряжения гауссовского шума) Я подозреваю, что делаю что-то не так, и хотел спросить, может ли кто-нибудь помочь с этим.
Вот код, объясняющий, что я пытаюсь сделать:
import numpy as np
from scipy.integrate import quad, simps
from matplotlib import pyplot as plt
def PDF(x, sigma=1, mu=0): # Gaussian normal distribution PDF
return 1/(np.sqrt(2*np.pi*sigma))*np.exp(-1/(2*sigma**2)*(x-mu)**2)
def PDFu(u, u_rms=1, u_mean=0):
return PDF(u, sigma=u_rms, mu=u_mean)
def PDFP(P):
return 2*PDFu(np.sqrt(P)) # substitute the input variable with the 'scaled' one
def probDensity(x, nbins): # calculate the probability density based on the input samples
distr, bins = np.histogram(x, nbins) # similar to plt.hist(density=True)
binWidth = bins[1]-bins[0]
binCenters = bins[:-1] binWidth/2
return distr/len(x)/binWidth, binCenters
npoints = 100000
rms = 1
u = np.random.normal(0, rms, npoints) # samples with Gaussian normal distribution
P = u**2 # square of the samples with Gaussian normal distribution - should follow chi-squared distribution?
nbins = 500
u_distr, u_bins = probDensity(u, nbins) # calculate PDF based on the samples
print('U_distr integral = ', simps(u_distr,u_bins)) # integrate the calculated PDF, should be 1
plt.plot(u_bins, u_distr)
us = np.linspace(-10, 10, 500)
PDFu_u = PDFu(us) # calculate the theoretical PDF
print('PDFu_u integral = ', quad(PDFu, -np.Inf, np.Inf)) # integral of the theoretical PDF, should be 1
plt.plot(us, PDFu_u)
nbins = 1000
P_distr, P_bins = probDensity(P, nbins) # calculate PDF based on the samples
print('P_distr integral = ', simps(P_distr, P_bins)) # integrate the calculated PDF, should be 1
plt.plot(P_bins, P_distr)
Ps = np.linspace(0, 8, npoints)
PDFP_P = PDFP(Ps) # calculate the theoretical PDF
plt.plot(Ps, PDFP_P)
print('PDFP_P integral = ', quad(PDFP, 0, np.Inf)) # integral of the theoretical PDF, should be 1
plt.show()
Теоретическое и моделируемое распределение вероятности нормальной случайной величины (u), похоже, хорошо совпадают, я использую это как проверку на вменяемость. Но разница существенна в случае переменной в квадрате, и я не могу понять, почему и как заставить их совпадать. Кстати, я пробовал различные правдоподобные коэффициенты масштабирования для теоретического распределения (например, 0,5, 2, sqrt (2)), но это не сработало, и я не понимаю, зачем мне это вообще нужно. Разве это не должно работать с простой заменой ‘P’ на ‘u’ в соответствии с формулой u = sqrt (P * R) [R = 1] и использованием нормального распределения ‘u’ для вычисления значения PDF для определенных ‘P’?
Я немного больше доверяю моделируемому распределению, и мне интересно, как правильно рассчитать теоретическое. Почему метод подстановки не работает?
Заранее благодарим вас за помощь!
Комментарии:
1. Я не совсем уверен, что вы сделали в программе, которую вы показали, но в любом случае X ^ 2 имеет распределение хи-квадрат на 1 д.е. Когда X является гауссовым со средним значением 0 и дисперсией 1.
Ответ №1:
Ваша теоретическая плотность для квадрата гауссова неверна. Вот calc. Если X является гауссовым, то для CDF $ F $ квадратной переменной $ Y = X ^ 2 $ имеем
$$ F(x) = P(Y<x) = P(X ^ 2 <x) = P(- sqrt{x} < X < sqrt{x}) = Phi( sqrt{x}) — Phi(- sqrt{x}) $$
где $ Phi $ — гауссовский CDF
итак, для PDF $ f (x) $ из $ Y $ мы дифференцируем это и получаем
$$ f(x) = F'(x) = (1/(2 sqrt{x})) Psi'(sqrt{x}) (1/(2sqrt{x})) Psi'(-sqrt{x}) = (1/(2 sqrt{x})) (psi(sqrt{x}) psi(- sqrt{x}) $$
где $ psi $ — это гауссовский PDF-файл
так что, по крайней мере, вам не хватает термина $(1/(2 sqrt{x})) $
Вот изображение формул, если это поможет 
Комментарии:
1. ой, извините, я думал, что здесь включен TeX. Есть ли что-то, что мне нужно сделать, чтобы показать формулы, кто-нибудь знает?
2. Спасибо! Объяснение имеет смысл. Коэффициент масштабирования 1 / (2 * sqrt (x)) между 2 PDF-файлами является ключевым отличием, которое необходимо для получения правильных теоретических значений. Я не понимаю (интуитивно), почему это необходимо, но математика работает.
Ответ №2:
Для справки, вот код с исправленным PDF, основанный на ответе Питербарга. Еще раз спасибо!
import numpy as np
from scipy.integrate import quad, simps
from matplotlib import pyplot as plt
def PDF(x, sigma=1, mu=0): # Gaussian normal distribution PDF
return 1/(np.sqrt(2*np.pi*sigma))*np.exp(-1/(2*sigma**2)*(x-mu)**2)
def PDFu(u, u_rms=1, u_mean=0):
return PDF(u, sigma=u_rms, mu=u_mean)
def PDFP(P):
return 1/(2*np.sqrt(P))*2*PDFu(np.sqrt(P)) # substitute the input variable with the 'scaled' one
def probDensity(x, nbins): # calculate the probability density based on the input samples
distr, bins = np.histogram(x, nbins) # similar to plt.hist(density=True)
binWidth = bins[1]-bins[0]
binCenters = bins[:-1] binWidth/2
return distr/len(x)/binWidth, binCenters
npoints = 100000
rms = 1
u = np.random.normal(0, rms, npoints) # samples with Gaussian normal distribution
P = u**2 # square of the samples with Gaussian normal distribution - should follow chi-squared distribution?
nbins = 500
u_distr, u_bins = probDensity(u, nbins) # calculate PDF based on the samples
print('U_distr integral = ', simps(u_distr,u_bins)) # integrate the calculated PDF, should be 1
plt.plot(u_bins, u_distr)
us = np.linspace(-10, 10, 500)
PDFu_u = PDFu(us) # calculate the theoretical PDF
print('PDFu_u integral = ', quad(PDFu, -np.Inf, np.Inf)) # integral of the theoretical PDF, should be 1
plt.plot(us, PDFu_u)
nbins = 1000
P_distr, P_bins = probDensity(P, nbins) # calculate PDF based on the samples
print('P_distr integral = ', simps(P_distr, P_bins)) # integrate the calculated PDF, should be 1
plt.plot(P_bins, P_distr)
Ps = np.linspace(0, 8, npoints)
PDFP_P = PDFP(Ps) # calculate the theoretical PDF
plt.plot(Ps, PDFP_P)
print('PDFP_P integral = ', quad(PDFP, 0, np.Inf)) # integral of the theoretical PDF, should be 1
plt.show()