7 Conceitos básicos sobre distribuições de probabilidade

O objetivo desta sessão é mostrar o uso de funções do R em problemas de probabilidade. Exercícios que podem (e devem!) ser resolvidos analiticamente são usados para ilustrar o uso do programa e alguns de seus recursos para análises numéricas.

Os problemas nesta sessão foram retirados do livro:
Bussab, W.O. & Morettin, P.A. Estatística Básica. $4^a$ edição. Atual Editora. 1987.

EXEMPLO 1 (adaptado de Bussab & Morettin, página 132, exercício 1) Dada a função

$$f(x) = \left\{ \begin{array}{ll} 2 \exp(-2x) & \mbox{ , se $x \geq 0$} \cr 0 & \mbox{ , se $x < 0$} \end{array} \right.$$

(a)

mostre que está função é uma f.d.p.

(b)

calcule a probabilidade de que $X > 1$

(c)

calcule a probabilidade de que $0.2 < X < 0.8$

Para ser f.d.p. a função não deve ter valores negativos e deve integrar 1 em seu domínio. Vamos começar definindo esta função como uma função no R para qual daremos o nome de $f1$. A seguir fazemos o gráfico da função. Como a função tem valores positivos para $x$ no intervalo de zero a infinito temos, na prática, para fazer o gráfico, que definir um limite em $x$ até onde vai o gráfico da função. Vamos achar este limite tentando vários valores, conforme mostram os comandos abaixo. O gráfico escolhido foi o produzido pelo comando plot(f1,0,5) e mostrado na Figura .

f1 <- function(x){
  fx <- ifelse(x < 0, 0, 2*exp(-2*x))
  return(fx)
}

plot(f1)
plot(f1,0,10)
plot(f1,0,5)

Figura : Gráfico da função de probabilidade do Exemplo 1.

Para verificar que a a integral da função é igual a 1 podemos usar a função integrate que efetua integração numérica. A função recebe como argumentos o objeto com a função a ser integrada e os limites de integração. Neste exemplo o objeto é f1 definido acima e o domínio da função é $[0, Inf]$. A saída da função mostra o valor da integral (1) e o erro máximo da aproximação numérica.

> integrate(f1, 0, Inf)
1 with absolute error < 5e-07

Para fazer cálculos pedidos nos itens (b) e (c) lembramos que a probabilidade é dada pela área sob a curva da função no intervalo pedido. Desta forma as soluções seriam dadas pelas expressões

$$\begin{eqnarray*} p_b & = & P(X > 1) = \int_1^\infty f(x) dx = \int_1^\infty 2\... ...) = \int_{0,2}^{0,8} f(x) dx = \int_{0.2}^{0.8} 2\,e^{-2x} dx \, \end{eqnarray*}$$

cuja representação gráfica é mostrada na Figura . Os comandos do R a seguir mostram como fazer o gráfico de função. O comando plot desenha o gráfico da função. Para destacar as áreas que correspondem às probabilidades pedidas vamos usar a função polygon. Esta função adiciona a um gráfico um polígono que é definido pelas coordenadas de seus vértices. Para sombrear a área usa-se o argumento density. Finalmente, para escrever um texto no gráfico usamos a função text com as coordenadas de posição do texto.

Figura : Probabilidades pedidas nos itens (b) e (c) do Exemplo 1.

> plot(f1,0,5)
> polygon(x=c(1,seq(1,5,l=20)), y=c(0,f1(seq(1,5,l=20))), density=10) 
> polygon(x=c(0.2,seq(0.2,0.8,l=20),0.8), y=c(0,f1(seq(0.2,0.8,l=20)), 0), col="gray") 
> text(c(1.2, 0.5), c(0.1, 0.2), c(expression(p[b],p[c])))

e como obter as probabilidades pedidas.

> integrate(f1, 1, Inf)
0.1353353 with absolute error < 2.1e-05
> integrate(f1, 0.2, 0.8)
0.4684235 with absolute error < 5.2e-15

EXEMPLO 2 (Bussab & Morettin, página 139, exercício 10) A demanda diária de arroz em um supermercado, em centenas de quilos, é uma v.a. $X$ com f.d.p.

$$f(x) = \left\{ \begin{array}{ll} \frac{2}{3}x \mbox{ , se $0 \leq... ... $ 1 \leq x < 3$} \cr 0 \mbox{ , se $x < 0$\ ou $x \geq 3$} \end{array} \right.$$ (1)

(a)

Calcular a probabilidade de que sejam vendidos mais que 150 kg.

(b)

Calcular a venda esperada em 30 dias.

(c)

Qual a quantidade que deve ser deixada à disposição para que não falte o produto em 95% dos dias?

Novamente começamos definindo um objeto do R que contém a função dada em .

Neste caso definimos um vetor do mesmo tamanho do argumento x para armazenar os valores de $f(x)$ e a seguir preenchemos os valores deste vetor para cada faixa de valor de $x$. A seguir verificamos que a integral da função é 1 e fazemos o seu gráfico mostrado na Figura .

> f2 <- function(x){
+  fx <- numeric(length(x))
+  fx[x < 0] <- 0
+  fx[x >= 0 & x < 1] <- 2*x[x >= 0 & x < 1]/3
+  fx[x >= 1 & x <= 3] <- (-x[x >= 1 & x <= 3]/3) + 1
+  fx[x > 3] <- 0
+  return(fx)
+ }

> integrate(f2, 0, 3) ## verificando que a integral vale 1
1 with absolute error < 1.1e-15

> plot(f2, -1, 4)     ## fazendo o gráfico da função

Figura : Gráfico da função densidade de probabilidade do Exemplo 2.

Agora vamos responder às questões levantadas. Na questão (a) pede-se a probabilidade de que sejam vendidos mais que 150 kg (1,5 centenas de quilos), portanto a probabilidade $P[X > 1,5]$. A probabilidade corresponde à área sob a função no intervalo pedido ou seja $P[X > 1,5] = \int_{1,5}^\infty f(x) dx$ e esta integral pode ser resolvida numericamente com o comando:

> integrate(f2, 1.5, Inf)
0.3749999 with absolute error < 3.5e-05

A venda esperada em trinta dias é 30 vezes o valor esperado de venda em um dia. Para calcular a esperança $E[X] = \int x f(x) dx$ definimos uma nova função e resolvemos a integral. A função integrate retorna uma lista onde um dos elementos ($value) é o valor da integral.

## calculando a esperança da variável
> ef2 <- function(x){ x * f2(x) }
> integrate(ef2, 0, 3)
1.333333 with absolute error < 7.3e-05

> 30 * integrate(ef2, 0, 3)$value ## venda esperada em 30 dias
[1] 40

Na questão (c) estamos em busca do quantil 95% da distribuição de probabilidades, ou seja o valor de $x$ que deixa 95% de massa de probabilidade abaixo dele. Este valor que vamos chamar de $k$ é dado por:

$$\int_0^k f(x) dx=0.95.$$

Para encontrar este valor vamos definir uma função que calcula a diferença (em valor absoluto) entre 0.95 e a probabilidade associada a um valor qualquer de $x$. O quantil será o valor que minimiza esta probabilidade. Este é portanto um problema de otimização numérica e para resolvê-lo vamos usar a função optimize do R, que recebe como argumentos a função a ser otimizada e o intervalo no qual deve procurar a solução. A resposta mostra o valor do quantil $x = 2.452278$ e a função objetivo com valor muito próximo de 0, que era o que desejávamos.

> f <- function(x) abs(0.95 - integrate(f2, 0, x)$value)  
> optimise(f, c(0,3))
$minimum
[1] 2.452278

$objective
[1] 7.573257e-08

A Figura  ilustra as soluções dos itens (a) e (c) e os comandos abaixo foram utilizados para obtenção destes gráficos.

par(mfrow=c(1,2), mar=c(3,3,0,0), mgp=c(2,1,0))
plot(f2, -1, 4)
polygon(x=c(1.5, 1.5, 3), y=c(0,f2(1.5),0), dens=10)

k <- optimise(f, c(0,3))$min
plot(f2, -1, 4)
polygon(x=c(0, 1, k, k), y=c(0,f2(1), f2(k), 0), dens=10)
text(c(1.5, k), c(0.2, 0), c("0.95", "k"), cex=2.5)

Figura : Gráficos indicando as soluções dos itens (a) e (c) do Exemplo 2.

Finalmente lembramos que os exemplos discutidos aqui são simples e não requerem soluções numéricas, devendo ser resolvidos analiticamente. Utilizamos estes exemplos somente para ilustrar a obtenção de soluções numéricas com o uso do R, que na prática deve ser utilizado em problemas mais complexos onde soluções analíticas não são triviais ou mesmo impossíveis.

7.1 Exercícios

1. (Bussab & Morettin, 5a edição, pag. 194, ex. 28)
   Em uma determinada localidade a distribuição de renda, em u.m. (unidade monetária) é uma variável aleatória $X$ com função de distribuição de probabilidade:
   
   $$f(x) = \left\{ \begin{array}{ll} \frac{1}{10}x + \frac{1}{10... ...} \cr 0 & \mbox{ se $x < 0$\ ou $x > 6$} \end{array} \right.$$
   
   
   1. mostre que $f(x)$ é uma f.d.p..
   2. calcule os quartis da distribuição.
   3. calcule a probabilidade de encontrar uma pessoa com renda acima de 4,5 u.m. e indique o resultado no gráfico da distribuição.
   4. qual a renda média nesta localidade?