1. Introdução

Considere o seguinte problema de caracterizar a distribuição (função cumulativa ou distribuição de probabilidade) de uma função de uma variável aleatória.

Se \(X\sim f_X(x)\), sendo \(f_X(x)\) (função de densidade de \(X\)) ou \(F_X(x)\) (função cumulativa de \(X\)) conhecidas, e \(Y=g(X)\), sendo \(g(x)\) uma função definida num domínio apropriado, como poderia caracterizar \(f_Y(y)\) e \(F_Y(y)\) ?

Há métodos algébricos disponíveis para resolver esse problema, como por exemplo o método da função cumulativa. Mas esses métodos, em muitos casos, não podem ser aplicados ou sua aplicação pode ser difícil.

Nesta nota técnica vamos resolver esse problema utilizando simulação Monte Carlo, uma técnica bem geral para esse tipo de situação, com resultados teoricamente justificados pela Lei dos Grandes Números.

2. Problema 1: encontrando \(\text{Pr}(Y\le 1)\)

Problema 1: Se \(X\sim\)Uniforme(0,2) e \(Y=X^2\), obtenha \(\text{Pr}(Y\le 1)\).

Para resolver por simulação vamos inicialmente definir \(n\) o número de simulações que desejamos considerar,

n<-10^7

Vamos agora definir um vetor \(x\) com os \(n\) valores simulados, considerando a distribuição Uniforme(0,2) definida para \(X\). A seguir vamos definir um vetor \(Y=X^2\). Essas operações seriam realizadas por

x<-runif(n,min=0,max=2)
y<-x^2

vamos “ver” os primeiros 5 valores de \(x\) e de \(y\):

x[1:5]
## [1] 0.01116459 0.65990361 0.32813745 0.60093848 0.35615793
y[1:5]
## [1] 0.000124648 0.435472772 0.107674184 0.361127055 0.126848473

Os valores que vai obter se repetir esse procedimento no seu computador vão ser diferentes, pois corresponderá a outro sorteio independente. Isso vale para todos os resultados apresentados nessa nota técnica: os resultados dependem dos valores que foram sorteados.

Podemos verificar o atendimento da condição \(Y\le 1\) para os primeiros 5 elementos do vetor \(y\), usando

y[1:5]<=1 # testando se cada elemento de y[1:5] é menor ou igual a 1
## [1] TRUE TRUE TRUE TRUE TRUE

que retorna um vetor com os resultados lógicos TRUE e FALSE do teste \(Y\le 1\).

Já poderiamos contar diretamente o número de resultados TRUE usando a função sum do R, que soma os elementos dentro de um vetor:

sum(y[1:5]<=1)
## [1] 5

Para efeito da soma, o TRUE tem valor 1 e o FALSE valor 0. E nesses primeiros 5 valores sorteados, teríamos uma frequência de simulações em que \(Y\le 1\) definida por

sum(y[1:5]<=1)/5
## [1] 1

Agora, para calcular a frequência de casos em que a condição \(Y\le 1\) é atendida em todas as \(n\) simulações, podemos usar:

ns<-sum(y<=1) # total se casos que atendem
ns/n          # frequência de atendimento -> probabilidade
## [1] 0.4998749

Podemos aferir a excelente qualidade da estimativa amostral, comparando esse o último valor acima com o valor exato (\(0{,}5\)) produzido pela caracterização algébrica da função cumulativa de \(Y\), que é disponível nesse caso: \[F_Y(y)=\frac{\sqrt{y}}{2}\ \text{e}\ \ F_Y(1)=\frac{1}{2}=0{,}5.\]

3. Problema 2: \(F_Y(y)\) como função do R

Problema 2: Se \(X\sim\)Uniforme(0,2) e \(Y=X^2\), defina uma função \(Fy(y,n)\) no R que emule em software a função cumulativa de \(Y\), por simulação, em que \(n\) define o número de simulações realizadas. Teste a função calculando com ela \(\text{Pr}(Y\le 1)=Fy(1)\), verificando se o valor obtido é próximo de \(0{,}50\).

Fy<-function(y,n=10^7){
 x<-runif(n,min=0,max=2)^2
 sum(x<=y)/n
}
Fy(1) # se não definir o n, usa o valor "default" definido na função
## [1] 0.5000811

4. Problema 3: Faça de gráfico de \(F_Y(y),\ \ y\in[-1,5]\)

Para fazer o gráfico, faremos inicialmente uma pequena alteração na definição da função Fy, para que ela permita avaliação vetorial, ou seja, ao entrar com um vetor de valores \(y\), ela já computa um vetor de probabilidades, retornando esse vetor como resultado. No código a seguir definimos a função e testamos com os valores no vetor c(1,2,3):

Fy<-function(y,n=10^7){
  z<-runif(n,min=0,max=2)^2
  x<-rep(0,length(y))
  for(i in 1:length(y))
     x[i]<-sum(z<=y[i])/n
  x
}

Fy(c(1,2,3))  # testando a função com o vetor c(1,2,3)
## [1] 0.4999790 0.7070584 0.8660449

Para desenhar uma função qualquer, que permita avaliação vetorial, os procedimentos gerais são simples. Primeiro estabeleça um número de pontos razoavelmente grande no domínio desejado, definido num vetor x (por exemplo). Segundo, avalie esse vetor na função, achando um vetor y (por exempo) com os resultados. Terceiro, “plote” os dois vetores, na ordem, com a função plot do R.

A seguir iremos apresentar o gráfico da função cumulativa usando recursos do comando plot. Caso não esteja familiarizado com esses procedimentos, recomendo dar uma olhada num tutorial que preparei sobre o uso da função plot para gráficos em geral, clicando neste link.

A seguir usaremos esses procedimentos recém apresentados para plotar Fy(y) num domínio [-1,5]

y<-seq(-1,5,0.1)
plot(y,Fy(y),type="l",col="blue")

Esses são os gráficos mais básicos que podemos fazer no R, usando o comando plot. Veja o link opções do comando plot, para mais detalhes sobre outras opções do comando plot que podem ser usadas.