Derivação de distribuições de probabilidade

Algumas vezes temos funções de densidade conhecidas para determinadas variáveis aleatórias (VA’s), mas queremos inferir qual seria o valor para uma outra VA. Por exemplo, conhecemos p_x(x) (densidade de probabilidade da variável ‘x’), sabemos que Y=2X (Y=g(X)) e precisamos encontrar p_y(y).

A abordagem mais intuitiva que conheço pra começar o problema é imaginar que temos uma função ‘g’, monótona crescente e com inversa g’(na falta de notação melhor), que mapeia elementos de X a Y (como qualquer função) e o segredo pra pegar a intuição do problema é lembrar que desta maneira, a probabilidade de qualquer intervalo [x0, x1] de p_x(x) é a mesma de qualquer intervalo [g(x0),g(x1)] em p_y(y). Deste modo podemos escrever que sempre a função acumulativa de Y, P_y(y) (que significa a probabilidade de que uma realização esteja no intervalo [-inf,y]) se relaciona com P_x(x) da seguinte maneira:

P_y(y)=P_x(g’(x))

Caso ‘g’ fosse monótona decrescente, estaríamos mapeando o intervalo ao contrário: [x0, x1] se relacionariam com [g(x1), g(x0)], a única mudança que isso acarreta é que teríamos P_y(y)=1-P_x(g’(x)).

Se quisermos calcular p_y(y), temos: p_y(y)=d P_y(Y)/dy=d P_x(g’(x))/dy

Exemplo numérico

Vamos fazer um exemplo numérico pra consolidar os conceitos. Imaginemos p_x(x) uma uniforme de 0 a 1 (U(0,1)) e Y=sqrt(X). Desta maneira temos a função densidade e a acumulativa de x nos gráficos seguintes:

Os gráficos foram gerados usando MATLAB com o código seguinte:

Bom e Y, como fica? Já vimos acima que P_y(y)=P_x(g’(y)) já que a raiz é uma função monótona crescente. Como estamos resolvendo tudo numericamente, o vetor com os valores de P_y(y) pode ser igual ao P_x(x), e usamos o intercâmbio entre as cumulativas onde P_y(g(x))=P_x(x). Assim pra plotar a CDF de Y, só temos que plotar a CDF de x, mas com o domínio em Y ou g(X).

Para obter a p_y(y) precisamos derivar a cumulativa, mas como estamos fazendo uma derivada numérica, temos que nos atentar ao domínio. O vetor ‘y’ contendo os possíveis valores da VA Y não está igualmente espaçado, portanto precisamos derivá-lo também. Neste caso temos alguns problemas numéricos com estas derivadas. Se calculamos usando a função ‘diff’ vamos estimar o que seria a derivada no meio dos pontos que utilizamos.

Exemplo: se aplicamos diff(x) vamos fazer um vetor [x(2)-x(1), x(3)-x(2),..] estas derivadas se aproximam mais do valor entre os pontos utilizados. Um outro problema é que temos um vetor menor que o inicial por uma unidade. Por outro lado, podemos usar a função ‘gradient’ que calcula a derivada usando 2 pontos, mas da seguinte maneira [x(2)-x(1),(x(3)-x(1))/2,…] e ao final temos um vetor de derivadas de mesmo tamanho que o de entrada. Nas extremidades ela calcula de uma forma e nos demais pontos, de outra. Você pode procurar por técnicas de derivação numérica e vai encontrar mais informações sobre isso.

O fato é que usando gradient(y) temos uma pequena descontinuidade no que seria a derivada correta, o que podemos observar no começo do gráfico de p_y(y).

Um material mais técnico sobre o tema pode ser encontrado aqui: https://www.engr.colostate.edu/~ramirez/ce_old/classes/ce522_ramirez/DerivedDistributions.pdf