PDF суммы гауссовых распределений с использованием БПФ

Question

Я пытаюсь вывести PDF из суммы независимых случайных величин. Сначала я хотел бы сделать это для простого случая: сумма гауссовских случайных величин.

Я был удивлен, увидев, что я не получаю гауссову функцию плотности, когда суммирую четное число гауссовских случайных величин. Я на самом деле получаю:

, который выглядит как две половины гауссовского распределения. С другой стороны, когда я суммирую нечетное число гауссовых распределений, я получаю правильное распределение:

ниже кода, который я использовал для получения результатов выше:

import numpy as np
from scipy.stats import norm
from scipy.fftpack import fft,ifft
import matplotlib.pyplot as plt
%matplotlib inline

a=10**(-15)
end=norm(0,1).ppf(a)
sample=np.linspace(end,-end,1000)
pdf=norm(0,1).pdf(sample)
plt.subplot(211)
plt.plot(np.real(ifft(fft(pdf)**2)))
plt.subplot(212)
plt.plot(np.real(ifft(fft(pdf)**3)))

Может ли кто-нибудь помочь мне понять, почему я получаю нечетные результаты для четных сумм гауссовых распределений?

Answer 1

Даже если ваш код создает гауссовское PDF с нулевым средним:

sample=np.linspace(end,-end,1000)
pdf=norm(0,1).pdf(sample)

БПФ не знает о sample и видит pdf только с выборками в 0, 1, 2, 3, ... 999. БПФ ожидает, что источником будет первая выборка сигнала. Для функции FFT ваш PDF не равен нулю, но имеет среднее значение 500.

Таким образом, здесь происходит добавление двух PDF-файлов со средним значением 500, что приводит к одному со средним значением 1000. А поскольку БПФ налагает периодичность на сигнал пространственной области, вы видите, что PDF выходит из графика справа и возвращается слева.

Добавление 3 PDF-файлов сдвигает среднее значение до 1500, что из-за периодичности равно 500, что означает, что оно заканчивается в том же месте, что и исходный PDF.

Решение состоит в том, чтобы сместить начало координат к первому образцу для БПФ и вернуть результат обратно:

from scipy.fftpack import fftshift, ifftshift
pdf2 = fftshift(ifft(fft(ifftshift(pdf))**2))

ifftshift смещает сигнал так, чтобы центральный сэмпл заканчивался на первом сэмпле, а fftshift смещает его обратно туда, где вы хотели его отобразить.

Но учтите, что при создании PDF исходная точка не является образцом, поэтому вышеприведенное не будет работать точно. Вместо этого используйте:

sample=np.linspace(end,-end,1001)
pdf=norm(0,1).pdf(sample)

При выборе 1001 выборки вместо 1000, ноль находится точно на средней выборке.

Answer 2

Используйте R!

library(ggplot2) 
f <- function(n) {
   x1 <- rnorm(n)
   x2 <- rnorm(n)
   X <- x1+x2
   return(ds)
}
ds.list <- lapply(10^(2:5),f)
ds <- Reduce(rbind,ds.list)
ggplot(ds,aes(X,fill = n)) + geom_density(alpha = 0.5) + xlab("")

Вот график распределения: