Qual é a distribuição de

17

Eu tenho quatro variáveis ​​independentes distribuídas uniformemente a,b,c,d , cada uma em [0,1] . Quero calcular a distribuição de (ad)2+4bc . Calculei a distribuição de u2=4bc para ser

f2(u2)=14lnu24
(daíu2(0,4]) e deu1=(ad)2para ser
f1(u1)=1u1u1.
Agora, a distribuição de uma somau1+u2é ( são também independentes) f u 1 + u 2 ( x ) = + - f 1 ( x - y ) f 2 ( y ) d y = - 1u1,u2porquey(0,4]. Aqui, ele deve serx>y,para que a integral seja igual afu1+u2(x)=-1
fu1+u2(x)=+f1(xy)f2(y)dy=14041xyxylny4dy,
y(0,4]x>yAgora eu o insiro no Mathematica e obtenho quefu1+u2(x)=1
fu1+u2(x)=140x1xyxylny4dy.
fu1+u2(x)=14[x+xlnx42x(2+lnx)].

Eu fiz quatro conjuntos independentes a,b,c,d consistindo em números cada e desenhei um histograma de ( a - d ) 2 + 4 b c :106(ad)2+4bc

insira a descrição da imagem aqui

e desenhou um gráfico de :fu1+u2(x)

insira a descrição da imagem aqui

Geralmente, o gráfico é semelhante ao histograma, mas no intervalo maioria é negativa (a raiz está em 2.27034). E o integrante da parte positiva é 0,77 .(0,5)0.77

Onde está o erro? Ou onde estou perdendo alguma coisa?

Edição: redimensionei o histograma para mostrar o PDF.

insira a descrição da imagem aqui

EDIT 2: Acho que sei onde está o problema no meu raciocínio - nos limites da integração. Como e x - y ( 0 , 1 ] , não posso simplesmente x 0. O gráfico mostra a região na qual tenho que integrar:y(0,4]xy(0,1]0x

insira a descrição da imagem aqui

Isso significa que eu tenho para y ( 0 , 1 ] (é por isso que parte do meu f estava correta), x x - 1 em y ( 1 , 4 ] e 40xy(0,1]fx1xy(1,4]x14y(4,5]

EDIÇÃO 3: Parece que o Mathematica PODE calcular as três últimas integrais com o seguinte código:

(1/4)*Integrate[((1-Sqrt[u1-u2])*Log[4/u2])/Sqrt[u1-u2],{u2,0,u1}, Assumptions ->0 <= u2 <= u1 && u1 > 0]

(1/4)*Integrate[((1-Sqrt[u1-u2])*Log[4/u2])/Sqrt[u1-u2],{u2,u1-1,u1}, Assumptions -> 1 <= u2 <= 3 && u1 > 0]

(1/4)*Integrate[((1-Sqrt[u1-u2])*Log[4/u2])/Sqrt[u1-u2],{u2,u1-1,4}, Assumptions -> 4 <= u2 <= 4 && u1 > 0]

que dá uma resposta correta :)

corey979
fonte
2
f1(u1)f2f1f2
f1f2f1f
Ao gerar esses histogramas para comparar com quantidades algébricas calculadas, dimensione o histograma para uma densidade válida (e sobreponha-os, se possível). Faça uma verificação semelhante para seus f1 e f2 para garantir que você tenha os direitos; se eles estão certos (ainda não vi nenhuma boa razão para suspeitar deles, mas é melhor checar novamente), o problema deve ser mais tarde.
Glen_b -Reinstala Monica

Respostas:

19

Muitas vezes, ajuda a usar funções de distribuição cumulativa.

Primeiro,

F(x)=Pr((uma-d)2x)=Pr(|uma-d|x)=1-(1-x)2=2x-x.

Próximo,

G(y)=Pr(4bcy)=Pr(bcy4)=0y/4dt+y/41ydt4t=y4(1log(y4)).

δ05(ad)2+4bcx=(ad)2Fy=4bcg=G

H(δ)=Pr((ad)2+4bcδ)=Pr(xδy)=04F(δy)g(y)dy.

Podemos esperar que isso seja desagradável - o PDF de distribuição uniforme é descontínuo e, portanto, deve produzir interrupções na definição deH- então é algo surpreendente que o Mathematica obtenha uma forma fechada (que não reproduzirei aqui). Diferenciando-o em relação aδdá a densidade desejada. É definido em partes dentro de três intervalos. Dentro0 0<δ<1,

H(δ)=h(δ)=18(8δ+δ(-(2+registro(16)))+2(δ-2δ)registro(δ)).

Dentro 1<δ<4,

h(δ)=14((δ+1)log(δ1)+δlog(δ)4δcoth1(δ)+3+log(4)).

And in 4<δ<5,

h(δ)=14(δ4δ4+(δ+1)log(4δ1)+4δtanh1((δ4)δδδδ4)1).

Figure

This figure overlays a plot of h on a histogram of 106 iid realizations of (ad)2+4bc. The two are almost indistinguishable, suggesting the correctness of the formula for h.


The following is a nearly mindless, brute-force Mathematica solution. It automates practically everything about the calculation. For instance, it will even compute the range of the resulting variable:

ClearAll[ a, b, c, d, ff, gg, hh, g, h, x, y, z, zMin, zMax, assumptions];
assumptions = 0 <= a <= 1 && 0 <= b <= 1 && 0 <= c <= 1 && 0 <= d <= 1; 
zMax = First@Maximize[{(a - d)^2 + 4 b c, assumptions}, {a, b, c, d}];
zMin = First@Minimize[{(a - d)^2 + 4 b c, assumptions}, {a, b, c, d}];

Here is all the integration and differentiation. (Be patient; computing H takes a couple of minutes.)

ff[x_] := Evaluate@FullSimplify@Integrate[Boole[(a - d)^2 <= x], {a, 0, 1}, {d, 0, 1}];
gg[y_] := Evaluate@FullSimplify@Integrate[Boole[4 b c <= y], {b, 0, 1}, {c, 0, 1}];
g[y_]  := Evaluate@FullSimplify@D[gg[y], y];
hh[z_] := Evaluate@FullSimplify@Integrate[ff[-y + z] g[y], {y, 0, 4}, 
          Assumptions -> zMin <= z <= zMax];
h[z_]  :=  Evaluate@FullSimplify@D[hh[z], z];

Finally, a simulation and comparison to the graph of h:

x = RandomReal[{0, 1}, {4, 10^6}];
x = (x[[1, All]] - x[[4, All]])^2 + 4 x[[2, All]] x[[3, All]];
Show[Histogram[x, {.1}, "PDF"], 
 Plot[h[z], {z, zMin, zMax}, Exclusions -> {1, 4}], 
 AxesLabel -> {"\[Delta]", "Density"}, BaseStyle -> Medium, 
 Ticks -> {{{0, "0"}, {1, "1"}, {4, "4"}, {5, "5"}}, Automatic}]
whuber
fonte
8
(+1), especially for reminding people that, instead say of density convolutions, "Often it helps to use cumulative distribution functions" -especially when they have such a simple form as here. And you were damn quick, also.
Alecos Papadopoulos
That looks like a neat solution that I'd love to accept - right after I understand it. I'm more a calculus man than a probabilist; at this moment I have three questions: i) how did you use the CDF to get F(x) and G(y), ii) why there's F and g under the integral for H, and iii) how do you from its form that the solution result will be piecewise?
corey979
(1) F and G are the CDFs. They are computed from the definition of a CDF, as indicated by the first equalities following their first appearances. The details should be apparent in the code I have inserted. (2) This is the convolution formula for a sum (more fully explained in a similar calculation at stats.stackexchange.com/a/144237). (3) I inserted a link to another thread about properties of uniform distributions.
whuber
7

Like the OP and whuber, I would use independence to break this up into simpler problems:

Let X=(ad)2. Then the pdf of X, say f(x) is:

enter image description here

Let Y=4bc. Then the pdf of Y, say g(y) is:

enter image description here

The problem reduces to now finding the pdf of X+Y. There may be many ways of doing this, but the simplest for me is to use a function called TransformSum from the current developmental version of mathStatica. Unfortunately, this is not available in a public release at the present time, but here is the input:

TransformSum[{f,g}, z]

which returns the pdf of Z=X+Y as the piecewise function:

enter image description here

Here is a plot of the pdf just derived, say h(z):

enter image description here

Quick Monte Carlo check

The following diagram compares an empirical Monte Carlo approximation of the pdf (squiggly blue) to the theoretical pdf derived above (red dashed). Looks fine.

enter image description here

wolfies
fonte