Cómo calcular el producto de Gauss variables en Matlab

Question

Quiero calcular la distribución de un producto de dos yo.yo.d. Gaussiano distribuido variables a y b. En principio, esto debería ser posible mediante la definición de una nueva variable x con una delta de dirac distribución

Para obtener la distribución de x (el producto de a y b), a y b tienen que ser marginados:

Desde la delta de dirac lazos de las variables a y b juntos a través de la variable x=ab, podemos saltar de uno integral y obtener

Hasta ahora tan bueno.

Si ahora se desea calcular esto en Matlab, tengo una enorme diferencia entre la solución con integral y un simple muestreo de la solución como aproximación.

Aquí está mi MWE:

mu_a = 2;
sig2_a = 1;

mu_b = 0;
sig2_b = 1;

res = 250;
xmin = -10;
xmax = 10;
x_range = linspace(xmin,xmax,res);

normal = @(x, mu, sig2) 1/sqrt(2*pi*sig2)*...
    exp((-1/2)*((x-mu)^2/sig2));

normal_a = @(a) arrayfun(normal,a,mu_a*ones(size(a)),sig2_a*ones(size(a)));
normal_b = @(b) arrayfun(normal,b,mu_b*ones(size(b)),sig2_b*ones(size(b)));

product_pdf1 = zeros(size(x_range));
product_pdf2 = zeros(size(x_range));
for i = 1:length(x_range)    
    x = x_range(i);

    normal_ab = @(a) normal_a(a).*normal_b(x./a);
    product_pdf1(i) = integral(normal_ab,-Inf,Inf);

    normal_ab = @(b) normal_a(x./b).*normal_b(b);
    product_pdf2(i) = integral(normal_ab,-Inf,Inf);
end

subplot(2,1,1,'replace'); hold all
plot(x_range, normal_a(x_range))
plot(x_range, normal_b(x_range))
plot(x_range, product_pdf1, '--')
plot(x_range, product_pdf2, '--')
axis tight
drawnow

% the sampling solution
N = 50000;
X_sequence = zeros(1,N);
for i=1:N
    a = mu_a + sqrt(sig2_a)*randn();
    b = mu_b + sqrt(sig2_b)*randn();
    x = a*b;
    X_sequence(i) = x;
end

subplot(2,1,2,'replace'); hold all
hist(X_sequence,250)
xlim([xmin xmax])

Tenga en cuenta que me calcula la integral de una sola vez y de una vez por b, el resultado es completamente diferente.

Dependiendo de los parámetros, uno, el otro o ambos de los cálculos son claramente erróneas. En mi ejemplo, el muestreo de la solución es mucho peakier que el mejor aspecto integral intento. El otro es completamente equivocado.

Ahora me pregunto: Es mi matemáticas mal? Es mi código de malo? O es una limitación de la de Matlab integral función (por ejemplo, tener problemas para la integración de más de discontinuidades)?

ACTUALIZACIÓN:

Tom descubrió un error en mis cálculos. He actualizado mi MWE en consecuencia y se normalizó el histograma para una mejor comparación:

mu_a = 2;
sig2_a = 1;

mu_b = 0;
sig2_b = 1;

res = 250;
xmin = -10;
xmax = 10;
x_range = linspace(xmin,xmax,res);

normal = @(x, mu, sig2) 1/sqrt(2*pi*sig2)*...
    exp((-1/2)*((x-mu)^2/sig2));

normal_a = @(a) arrayfun(normal,a,mu_a*ones(size(a)),sig2_a*ones(size(a)));
normal_b = @(b) arrayfun(normal,b,mu_b*ones(size(b)),sig2_b*ones(size(b)));

product_pdf1 = zeros(size(x_range));
product_pdf2 = zeros(size(x_range));
for i = 1:length(x_range)    
    x = x_range(i);

    normal_ab = @(a) (1./a).*normal_a(a).*normal_b(x./a);
    product_pdf1(i) = integral(normal_ab,-Inf,Inf);

    normal_ab = @(b) (1./b).*normal_a(x./b).*normal_b(b);
    product_pdf2(i) = integral(normal_ab,-Inf,Inf);
end

subplot(2,1,1,'replace'); hold all
plot(x_range, normal_a(x_range))
plot(x_range, normal_b(x_range))
plot(x_range, product_pdf1, '--')
plot(x_range, product_pdf2, '--')
axis tight
drawnow

% the sampling solution
N = 50000;
X_sequence = zeros(1,N);
for i=1:N
    a = mu_a + sqrt(sig2_a)*randn();
    b = mu_b + sqrt(sig2_b)*randn();
    x = a*b;
    X_sequence(i) = x;
end

subplot(2,1,2,'replace'); hold all
hres = 250; % histogram resolution
[B, X] = hist(X_sequence,hres);
dx = diff(X(1:2));
plot(x_range, product_pdf1,'--')
plot(X, B/sum(B*dx),'-.');
plot(x_range, besselk(0,abs(x_range))/pi)
xlim([xmin xmax])

Pero, extrañamente, las soluciones aún no coinciden: Uno de los integral soluciones ahora está completamente apagado, pero la otra (parte inferior de la parcela, de punto-rojo discontinuo) ahora está en realidad más cerca de la toma de muestras de la solución (discontinua verde). El BesselK no coincide, es mucho peakier.

Más sugerencias?

Answer 1

Usted cometió un error en la derivación. Vamos a empezar a partir de esta integral: $$ p(x) = \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \delta(ab-x) p(a) p(b) da ~ db $$ Para quitar el $\delta$, usted necesita para hacer un cambio de variable de$b$$c=ab$, por lo que $dc = a ~ db$. La integral se convierte en: $$ p(x) = \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \delta(c-x) p(a) p(c/a) \frac{1}{|a|} da ~ dc \\ = \int_{-\infty}^{\infty} p(a) p(x/a) \frac{1}{|a|} da $$ Al $p(a)$ es normal estándar, $p(x) = \frac{1}{\pi} {\rm BesselK}(0,|x|)$. Usted puede calcular esto en Matlab con el código besselk(0,abs(x))/pi.

Answer 2

Su código de MATLAB se ve bien. También está de acuerdo con mi propia implementación de las ecuaciones:

mu2 = 2;
[a, x] = ndgrid(linspace(-10, 10, 250));
clf; hold all
hist(randn(1,1e6) .* (randn(1, 1e6) + mu2), x(1,:));
y = sum(exp(-0.5.*(a-mu2).^2).*exp(-0.5.*(x./a).^2))';
y = y * (1e6 / sum(y));
plot(x(1,:)', y, 'r-');
y = sum(exp(-0.5.*a.^2).*exp(-0.5.*((x./a)-mu2).^2))';
y = y * (1e6 / sum(y));
plot(x(1,:)', y, 'g-');
xlim(x(1,[1 end]));

El hecho de que nuestros dos diferentes enfoques para calcular las integrales de acuerdo, y que las dos integrales de acuerdo tanto el uno con el otro y con la distribución empírica, me sugiere que el problema está en la matemática de la derivación. Sin embargo, yo no podía ver, pero que la pregunta no es apropiado para este intercambio de todos modos. Te sugiero que publiques en el de las Matemáticas o de la Cruz Validado intercambios.

Answer 3

Creo que su derivación está en mal estado. Echa un vistazo a este Wiki.

Aquí la implementación en MATLAB:

 f1=@(x,z)normpdf(x).*normpdf(z./x)./abs(x)
 f=@(z)integral(@(x)f1(x,z),-10,10)
 >> f(1)
 ans =
 0.1340

He aquí cómo usted puede simular:

 samp=prod(randn(2,100000));
 hist(samp,100);

Distribución del uso del accesorio de la herramienta de ajuste no paramétrico de la distribución de comparar los teóricos de la distribución:

 dfittool(samp')