require(gtools)
require(tidyverse)
require(ggforce)
options(kableExtra.latex.load_packages = TRUE)
require(kableExtra)
require(mosaicCalc)
require(gridExtra)
require(numDeriv)
require(DT)
require(pracma)
#Para numerar tabelas e figuras:
require(captioner)
fig_nums <- captioner(prefix = "Figura")
table_nums <- captioner(prefix = "Tabela")
quadro_nums<- captioner(prefix="Quadro")
Com frequência, conhecemos a distribuição de probabilidade de uma variável aleatória e estamos interessados em determinar a distribuição de alguma função dessa variável. Por exemplo, suponha que conheçamos a distribuição de \(X\) e queiramos obter a distribuição de \(Y=g(X)\). A densidade de \(Y\) pode ser calculada por derivação de \(F_Y(y)\): \[\small F_Y(y) = \int\limits_{-\infty}^{y}g(w)dw.\] Nesse caso, \(g\) é a densidade da variável \(X\).
Exemplo 1: Seja \(X\) uma variável aleatória contínua com densidade \(f(x)=\frac{1}{2}I_{\left[1,3\right]}(x)\). Encontre a função distribuição e a função densidade de \(Z=e^{X}\).
Solução: Vemos que a variável \(e^X\) assume valores de \(e \approx 2.7183\) até \(e^3\approx 20.0855\):
x = seq(1,3,0.5)
z = exp(x)
data=data.frame(x,z) %>%
`colnames<-`(c("x","z=exp(x)"))%>%
print()
## x z=exp(x)
## 1 1.0 2.718282
## 2 1.5 4.481689
## 3 2.0 7.389056
## 4 2.5 12.182494
## 5 3.0 20.085537
E a função de densidade acumulada de X: \[ \small \begin{array}{lll} F_X(x)=\left\{ \begin{array}{lll} 0 \mbox{, se x < 1}\\ \int\limits_1^x \frac{1}{2}du = \left.\frac{u}{2}\right|_{1}^{x} = \frac{x-1}{2} \mbox{, se } 1 \leq x \leq 3\\ 1 \mbox{, se } x >3 \end{array} \right. \end{array} \small \begin{array}{lll} \mbox{Pelo método da acumulada:}\\ P(Z<z)=P(e^X<z)=P(X<\log(z)) \end{array} \] \[ \small \begin{array}{lll} &&F_Z(z) = \left\{ \begin{array}{lll} 0, &\mbox{ se log(z) < 1, ou seja, z < e}\\ \frac{\log(z)-1}{2} &\mbox{ se } 1 \leq \log(z) \leq 3, \mbox{ ou seja}, e \leq z \leq e^3\\ 1, &\mbox{ caso contrário.} \end{array} \right.\\ &&\mbox{Obtendo a densidade de Y:}\\ &&f_Z(z)=\frac{\partial F_Z(z)}{\partial z}=\frac{1}{2z}I_{[e,e^3]}(z) \\ &&\mbox{ não conhecida no Apêndice Mood} \end{array} \]
Gráficos no R:
x = seq(0,4,0.1)
fx = ifelse(x<1,0,
ifelse(x<3,1/2,0))
Fx = ifelse(x<1,0,
ifelse(x<3,1/2*(x-1),1))
data = data.frame(x=x[-length(x)],fx=fx[-length(x)], xend=x[-1], fxend=fx[-length(x)],Fx=Fx[-length(x)])
a=ggplot(data, aes( x=x, y=fx ,xend=xend,yend=fxend)) +
geom_segment(size=1.2,color="orange")+
geom_point( aes( x=1, y=0 ), shape=1, size=3,color="orange")+
geom_point( aes( x=3, y=0 ), shape=1, size=3,color="orange") +
labs(x="x",y="f(x)",title="densidade de X ~ U(0,3)")
b=ggplot(data, aes( x=x, y=Fx)) +
geom_line(size=1.2,color="orange")+
labs(x="x",y="F(x)",title="densidade acumulada de X ~ U(0,3)")
grid.arrange(a, b, ncol = 2, nrow = 1)
Figura 1: Gráficos da função de densidade e densidade acumulada de X
De fato, a área sob a curva é igual a 1:
Area=antiD(1/(2*z) ~ z)
Area(exp(3))-Area(exp(1))
## [1] 1
Gráficos no R:
z = seq(0,23,0.1)
fz = ifelse(z<exp(1),0,
ifelse(z<exp(3),1/(2*z),0))
Fz = ifelse(z<exp(1),0,
ifelse(z<exp(3),(log(z)-1)/2,1))
data = data.frame(z=z[-length(z)],fz=fz[-length(z)], zend=z[-1], fzend=fz[-length(z)],Fz=Fz[-length(z)])
a=ggplot(data, aes( x=z, y=fz ,xend=zend,yend=fzend)) +
geom_segment(size=1.2,color="orange")+
geom_point( aes( x=exp(1), y=0 ), shape=1, size=3,color="orange")+
geom_point( aes( x=exp(3), y=0 ), shape=1, size=3,color="orange") +
labs(x="x",y="f(x)",title="densidade de Z")
b=ggplot(data, aes( x=z, y=Fz)) +
geom_line(size=1.3,color="orange")+
labs(x="z",y="F(z)",title="densidade acumulada de Z")
grid.arrange(a, b, ncol = 2, nrow = 1)
Figura 2: Gráficos das funções de densidade e densidade acumulada de Z
Exemplo 2: Seja \(X\) uma variável aleatória contínua com densidade \(f\). Obtenha a densidade de uma variável aleatória \(Z=X^2\). Pelo método da acumulada: \[\small \begin{array}{lll} P(Z<z)=P(X^2<z)=P(-\sqrt{z}<X<\sqrt{z})=P(X<\sqrt{z})-P(X<-\sqrt{z})\\ \Rightarrow F_Z(z) = F_X(\sqrt{z})-F_X(-\sqrt{z})\\ f_Z(z)=\frac{\partial F_Z(z)}{\partial z}=\frac{\partial F_X(\sqrt{z})}{\partial z}-\frac{\partial F_X(-\sqrt{z})}{\partial z}\\ =f_X(\sqrt{z}).\frac{\partial (\sqrt{z})}{\partial z}-f_X(-\sqrt{z}).\frac{\partial (-\sqrt{z})}{\partial z}\\ = f_X(\sqrt{z}).\frac{1}{2\sqrt{z}}+f_X(-\sqrt{z}).\frac{1}{2\sqrt{z}}\\ = \frac{1}{2\sqrt{z}} \left[f_X(\sqrt{z})+f_X(-\sqrt{z})\right] \end{array} \]
Exemplo prático considerando \(X \sim \mbox{Gamma}(2,4)\) A área sob a curva de Z é igual a 1:
a=2
b=4
f=function(z){
fx1=dgamma(sqrt(z),shape=a,scale=b)
fx2=dgamma(-sqrt(z),shape=a,scale=b)
1/(2*sqrt(z))*(fx1+fx2)
}
Area=integrate(f, lower = 0, upper = Inf)
Area
## 1 with absolute error < 4.5e-05
a=2
b=4
x = seq(0.01,25,0.1)
z = x^2
fx = dgamma(x,shape=a,scale=b)
fx1=dgamma(sqrt(z),shape=a,scale=b)
fx2=dgamma(-sqrt(z),shape=a,scale=b)
fz = 1/(2*sqrt(z))*(fx1+fx2)
dados=data.frame(x,z,fx,fz)
a=ggplot(dados, aes( x=x, y=fx)) +
geom_line(size=1.2,color="orange") +
labs(x="x",y="f(x)",title="densidade de X ~ Gamma(2,4)")
b=ggplot(dados, aes( x=z, y=fz)) +
geom_line(size=1.2,color="orange") +
labs(x="x",y="f(x)",title="densidade de Z=X^2")
grid.arrange(a, b, ncol = 2, nrow = 1)
Figura 3: Gráficos das funções de densidade de X e densidade Z
Exemplo prático considerando \(X \sim U(-3,-1)\) A área sob a curva de Z é igual a 1: considere integração numérica pela praticidade!
a=-3
b=-1
f=function(z){
fz1=dunif(sqrt(z),min=a,max=b)
fz2=dunif(-sqrt(z),min=a,max=b)
1/(2*sqrt(z))*(fz1+fz2)
}
Area=integrate(f, lower = 0, upper = Inf)
Area
## 1 with absolute error < 1.2e-06
a=-3
b=-1
x = seq(-4,0,0.05)
z = seq(0,10,0.05)
fx = dunif(x,min=a,max=b)
fz1=dunif(sqrt(z),min=a,max=b)
fz2=dunif(-sqrt(z),min=a,max=b)
fz = 1/(2*sqrt(z))*(fz1+fz2)
data1 = data.frame(x=x[-length(x)],fx=fx[-length(x)], xend=x[-1], fxend=fx[-length(x)])
data2 = data.frame(z=z[-length(z)],fz=fz[-length(z)], zend=z[-1], fzend=fz[-length(z)])
a=ggplot(data1, aes( x=x, y=fx,xend=xend,yend=fxend)) +
geom_segment(size=1.2,color="orange")+
geom_point( aes( x=-3, y=0 ), shape=1, size=3,color="orange")+
geom_point( aes( x=-1, y=0), shape=1, size=3,color="orange") +
labs(x="x",y="f(x)",title="densidade de X ~ U(-3,-1)")
b=ggplot(data2, aes( x=z, y=fz,xend=zend,yend=fzend)) +
geom_segment(size=1.2,color="orange")+
geom_point( aes( x=1, y=0 ), shape=1, size=3,color="orange")+
geom_point( aes( x=9, y=0 ), shape=1, size=3,color="orange") +
labs(x="z",y="f(z)",title="densidade de Z=X^2")
grid.arrange(a, b, ncol = 2, nrow = 1)
## Warning: Removed 1 rows containing missing values (geom_segment).
Figura 4: Gráficos das funções de densidade de X e densidade Z
Exemplo 3: Seja \(W\) uma variável aleatória contínua com densidade Normal \(N(0,\sigma^2)\). Obtenha a densidade da variável aleatória \(W=X^2\).
Solução: Vemos que a variável \(X^2\) assume valores positivos:
Gráficos no R para \(N(0,4)\):
media=0
desvio=2
x = seq(-10,10,0.1)
fx = dnorm(x,mean=media,sd=desvio)
Fx = pnorm(x,mean=media,sd=desvio)
dados=data.frame(x,fx,Fx)
a=ggplot(dados, aes( x=x, y=fx)) +
geom_line(size=1.2,color="orange") +
labs(x="x",y="f(x)",title="densidade de X ~ N(0,4)")
b=ggplot(dados, aes( x=x, y=Fx)) +
geom_line(size=1.2,color="orange") +
labs(x="x",y="f(x)",title="densidade acumulada de X ~ N(0,4)")
grid.arrange(a, b, ncol = 2, nrow = 1)
Figura 5: Gráficos das funções de densidade e densidade acumulada de X
Pelo método da acumulada: \[\small P(W<w)=P(X^2<w)=P(-\sqrt{w}<X<\sqrt{w})=P(X<\sqrt{w})-P(X<\sqrt{w})\\ \Rightarrow F_W(w) = F_X(\sqrt{w})-F_X(-\sqrt{w})\\ \mbox{ lembrando que a acumulada da normal não tem forma fechada! } \]
Mesmo assim, podemos obter a densidade de W - e verificamos que W não pode assumir valores nulos: \[\small \begin{array}{ll} f_W(w)&=\frac{\partial F_W(w)}{\partial w}\\ &=\frac{\partial F_X(\sqrt{w})}{\partial w}-\frac{\partial F_X(-\sqrt{w})}{\partial w}\\ &=f_X(\sqrt{w}).\frac{\partial (\sqrt{w})}{\partial w}-f_X(-\sqrt{w}).\frac{\partial (-\sqrt{w})}{\partial w} \mbox{ e com a simetria da normal: }\\ &=f_X(\sqrt{w}).\left(\frac{1}{2\sqrt{w}}+ \frac{1}{2\sqrt{w}} \right)\\ &=\frac{1}{\sqrt{2\pi} \sigma \sqrt{w}}.\exp\left[-\frac{1}{2\sigma^2}w\right]\\ &\mbox{ não conhecida no apêndice Mood, exceto para }\sigma^2=1 \\ &\mbox{ onde temos a distribuição} X^2 \mbox{ com 1 grau de liberdade!} \end{array} \] De fato, a área sob a curva é igual a 1:
sigma=2
f=function(w){
1/(sqrt(2*pi*w)*sigma)*exp(-1/(2*sigma^2)*w)
}
Area=integrate(f, lower = 0, upper = Inf)
Area
## 1 with absolute error < 5.2e-05
Gráficos no R considerando \(\sigma^2=4\):
media=0
desvio=2
w = seq(0.01,10,0.1)
fw = 1/(sqrt(2*pi*w)*desvio)*exp(-1/(2*desvio^2)*w)
Fw = pnorm(sqrt(w),mean=media,sd=desvio)-pnorm(-sqrt(w),mean=media,sd=desvio)
dados=data.frame(w,fw,Fw)
a=ggplot(dados, aes( x=w, y=fw)) +
geom_line(size=1.2,color="orange") +
labs(x="w",y="f(w)",title="densidade de W")
b=ggplot(dados, aes( x=w, y=Fw)) +
geom_line(size=1.2,color="orange") +
labs(x="w",y="F(w)",title="densidade acumulada de W")
grid.arrange(a, b, ncol = 2, nrow = 1)
Figura 6: Gráficos das funções de densidade e densidade acumulada de W
Teorema: Seja \(\varphi\) uma função diferenciável, estritamente crescente ou estritamente decrescente em um intervalo \(I\), e sejam \(\varphi(I)\) o contradomínio de \(\varphi\) e \(\varphi^{-1}\) a função inversa de \(\varphi\). Seja \(X\) uma variável aleatória contínua de densidade \(f\) tal que \(f(x) = 0\) para \(x\notin I\). Então \(Y = \varphi(X)\) tem densidade \(g\) dada por \(g(y)=0\) para \(y\notin \varphi(I)\) e \[\small g(y) = f(\varphi^{-1}(y))\left|\frac{d}{dy}\varphi^{-1}(y)\right|I_{\mathscr{D}}(y), \mathscr{D} \quad \mbox{é o conjunto domínio de Y}. \]
ou de forma equivalente: \[\small g(y) = f(x)\left|\frac{dx}{dy}\right|I_{\mathscr{D}}(y), \mathscr{D} \quad \mbox{é o conjunto domínio de Y} \quad \mbox{e} \quad x = \varphi^{-1}(y). \]Exemplo 4: Seja \(X\) uma variável aleatória tendo uma densidade exponencial de parâmetro \(\lambda\). Obtenha a densidade de \(Z = X^{1/\beta},\) onde \(\beta > 0.\)
Solução: Vemos que a variável \(Z\) assume somente valores positivos:
x = seq(0.1,2,0.3)
beta1=2
beta2=1/2
z1 = x^(1/beta1)
z2 = x^(1/beta2)
data.frame(x,z1,z2) %>%
`colnames<-`(c("x","z1=x^{1/beta} sendo beta=2","z2=x^{1/beta} sendo beta =1/2"))%>%
datatable(cap=quadro_nums("quadro1","Alguns valores que X e Z assumem"),
options = list(columnDefs=list(list(class="dt-center",targets=list(1,2,3)),
list(width = '20%', targets = list(1,2,3)))))
Pelo método do Jacobiano: \[\small \begin{array}{llll} & w=x^{\frac{1}{\beta}} \Rightarrow x=w^\beta \Rightarrow \left|\frac{dx}{dw}\right|=\beta w^{\beta-1}\\ f_W(w)& =f_X(x)\left|\frac{dx}{dw}\right|.I_{(0,\infty)}(w) \mbox{ , ou seja,}\\ & =f_X(w^\beta)\beta w^{\beta-1}.I_{(0,\infty)}(w) \\ & \mbox{ lembre-se que a função é de w então tiramos os termos x da fórmula!}\\ & = \lambda \beta w^{\beta-1} \exp(-\lambda w^{\beta}) .I_{(0,\infty)}(w) \\ &\mbox{ que corresponde a Weibull no apêndice Mood!}\\ &\therefore W \sim \mbox{ Weibull} (\lambda,\beta) \end{array} \]
Gráficos no R considerando \(\lambda=2\) e \(\beta=2\):
lambda=2
beta=2
x= seq(0.001,2,0.01)
w= seq(0.001,2,0.01)
fx = dexp(x,rate=lambda)
fw = dweibull(w,scale=1/lambda,shape=beta)
dados=data.frame(x,w,fx,fw)
a=ggplot(dados, aes( x=x, y=fx)) +
geom_line(size=1.2,color="orange") +
labs(x="x",y="f(x)",title="densidade de X: Exp(2)")
b=ggplot(dados, aes( x=w, y=fw)) +
geom_line(size=1.2,color="orange") +
labs(x="w",y="fw",title="densidade de W: Weibull(2,2)")
grid.arrange(a, b, ncol = 2, nrow = 1)
Figura 7: Gráficos das funções de densidade de X e densidade de W
Sejam \(X_1\) e \(X_2\) variáveis aleatórias contínuas com função densidade de probabilidade conjunta \(f_{X_1,X_2}(x_1,x_2)\). Além disso, sejam \(Y_1=g_1(x_1,x_2)\) e \(Y_2=g_2(x_1,x_2)\), tais que \(g_1\) e \(g_2\) satisfaçam as seguintes condições:
Nessas condições, pode-se mostrar que as variáveis aleatórias \(Y_1\) e \(Y_2\) são conjuntamente contínuas com função densidade conjunta dada por \[\small f_{Y_1,Y_2}(y_1,y_2) = f_{X_1,X_2}(x_1,x_2)\left|J(y_1,y_2)\right|I_{\mathscr{D}}(y_1,y_2), \] em que \(x_1 = g_1^{-1}(y_1,y_2)\), \(x_2 = g_2^{-1}(y_1,y_2)\) e \(\mathscr{D}\) é o domínio de \((Y_1,Y_2)\).
Exemplo 5: Sejam \(X\) e \(Y\) variáveis aleatórias conjuntamente contínuas com função densidade de probabilidade \(f_{X,Y}\). Sejam \(Z = X+Y\) e \(W = X-Y\).
Solução: \[ \small \begin{array}{ll} \mbox{Sejam } \left\{ \begin{array}{ll} z =g(x,y) = x+y \\ w = h(x,y) = x-y \end{array} \right. \Rightarrow \left\{ \begin{array}{ll} z+w=2x \rightarrow x=\frac{z+w}{2}\\ z-w=2y \rightarrow y=\frac{z-w}{2} \end{array} \right. \end{array}\\ \begin{array}{ll} \mbox{ e o det. do Jacobiano:}\\ J(z,w) = \left| \begin{array}{cc} \frac{dx}{dz} & \frac{dx}{dw} \\ \frac{dy}{dz} & \frac{dy}{dw} \end{array} \right| \end{array} =\left| \begin{array}{ll} \frac{1}{2} & \frac{1}{2} \\ \frac{1}{2} & -\frac{1}{2} \end{array} \right| =\left|-\frac{1}{4}-\frac{1}{4}\right|=\frac{1}{2} \] \[ \small \begin{array}{lll} \Rightarrow f_{Z,W}(z,w) &&= J(z,w).f_{X,Y}(x,y)I_{\mathscr{D}}(z,w)\\ &&=\frac{1}{2}.f_{X,Y}(\frac{z+w}{2},\frac{z-w}{2}).I_{\mathscr{D}}(z,w), \end{array} \] onde \(I_{\mathscr{D}}(z,w)\) é o indicador do domínio, ou seja, os valores que z e w podem assumir.
Solução: \[ \small f_{Z,W}(z,w) = \left\{ \begin{array} {ll} \frac{\lambda_1\lambda_2}{2}\exp\left\{-\lambda_1\left(\frac{z+w}{2}\right) -\lambda_2\left(\frac{z-w}{2}\right)\right\}, \mbox{ se } 0<\frac{z+w}{2}<\infty \mbox{ e } 0<\frac{z-w}{2}<\infty;\\ 0, \mbox{ caso contrário.} \end{array} \right.\\ \]
Quando a função densidade conjunta das \(n\) variáveis aleatórias \(X_1, X_2, \ldots, X_n\) é dada e queremos calcular a função densidade conjunta de \(Y_1, Y_2, \ldots, Y_n\), onde \[\small Y_1 = g_1(X_1,\ldots,X_n), \quad Y_2 = g_2(X_1,\ldots,X_n), \ldots, Y_n = g_n(X_1,\ldots,X_n),\] a abordagem é a mesma. Supomos que as equações \(y_1 = g_1(x_1,x_2,\ldots,x_n), y_2 = g_2(x_1,x_2,\ldots,x_n), \ldots, y_n = g_n(x_1,x_2,\ldots,x_n)\) tenham soluções únicas \(x_1 = g_1^{-1}(y_1,y_2,\ldots,y_n), \ldots, x_n = g_n^{-1}(y_1,y_2,\ldots,y_n)\) cuja função densidade conjunta das variáveis aleatórias \(Y_i\) é dada por \[\small f_{Y_1,\ldots,Y_n}(y_1,\ldots,y_n) = f_{X_1,\ldots,X_n}(x_1,\ldots,x_n)\left|J(y_1,\ldots,y_n)\right|I_{\mathscr{D}}(y_1,\ldots,y_n), \] em que \(x_i = g_i^{-1}(y_1,\ldots,y_n), i = 1,2,\ldots,n\) e
\[\small J(y_1,\ldots,y_n) = \left| \begin{array}{llll} \frac{\partial x_1}{\partial y_1} \frac{\partial x_1}{\partial y_2}\ldots \frac{\partial x_1}{\partial y_n}\\ \frac{\partial x_2}{\partial y_1} \frac{\partial x_2}{\partial y_2}\ldots \frac{\partial x_2}{\partial y_n}\\ \vdots \\ \vdots \ldots\\ \frac{\partial x_n}{\partial y_1} \frac{\partial x_n}{\partial y_2}\ldots \frac{\partial x_n}{\partial y_n} \end{array} \right| \neq 0. \]
Exemplo 6: Considere o vetor aleatório \(\textbf{X}\) com densidade \[\small f_\textbf{X}(\textbf{x}) = \lambda^3\mbox{e}^{-\lambda(x_1+x_2+x_3)}I_{(0,\infty)}(x_1)I_{(0,\infty)}(x_2)I_{(0,\infty)}(x_3), \lambda>0.\] Qual a densidade conjunta de \(\textbf{Y} = (Y_1,Y_2,Y_3)\) tal que \[\small \begin{array}{lll} Y_1 &= X_1;\\ Y_2 &= X_1+X_2;\\ Y_3 &= X_1+X_2+X_3. \end{array} \]
Solução: \[ \small \begin{array}{ll} \mbox{Sejam } \left\{ \begin{array}{ll} y_1 = x_1 \\ y_2 = x_1 + x_2\\ y_3 = x_1 + x_2 + x_3 \end{array} \right. \Rightarrow \left\{ \begin{array}{ll} x_1=y_1\\ x_2=y_2-y_1\\ x_3=y_3-y_2 \end{array} \right. \end{array}\\ \begin{array}{ll} \mbox{ e o det. do Jacobiano:}\\ J(y_1,y_2,y_3) = \left| \begin{array}{llll} \frac{\partial x_1}{\partial y_1} \frac{\partial x_1}{\partial y_2}\frac{\partial x_1}{\partial y_3}\\ \frac{\partial x_2}{\partial y_1} \frac{\partial x_2}{\partial y_2} \frac{\partial x_2}{\partial y_3}\\ \frac{\partial x_n}{\partial y_1} \frac{\partial x_n}{\partial y_2} \frac{\partial x_n}{\partial y_3} \end{array} \right| =\left| \begin{array}{ccc} 1 & 0 & 0 \\ -1 & 1 & 0\\ 0 & -1 & 1 \end{array} \right| =1 \end{array}\\ \begin{array}{lll} \Rightarrow f_{Y_1,Y_2,Y_3}(y_1,y_2,y_3)&=J(y_1,y_2,y_3).f_{X_1,X_2,X_3}(x_1,x_2,x_3).I_{\mathscr{D}}(y_1,y_2,y_3)\\ &=\lambda^3\mbox{e}^{-\lambda(y_3)}I_{(0,\infty)}(y_1)I_{(y_1,\infty)}(y_2)I_{(y_2,\infty)}(y_3), \lambda>0. \end{array} \]
O determinante da matriz jacobiana é igual a 1:
J=matrix(c(1,0,0,-1,1,0,0,-1,1),ncol=3,nrow=3,byrow=T)
J
## [,1] [,2] [,3]
## [1,] 1 0 0
## [2,] -1 1 0
## [3,] 0 -1 1
det(J)
## [1] 1