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.
Note que há uma edição mais nova: (\(5^a\) edição, 2003 - Ed. Saraiva)
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 . \]
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 e mostrado na Figura 16 foi o produzido pelo comando plot(f1,0,5).
> f1 <- function(x){ + fx <- ifelse(x < 0, 0, 2*exp(-2*x)) + return(fx) + } > plot(f1)
> plot(f1,0,10)
> plot(f1,0,5)
Para verificar que 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, \infty ]\). 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
cuja representação gráfica é mostrada na Figura 17. 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.
> 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 para obter as probabilidades pedidas usamos integrate().
> 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.
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\).
> 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) + }
A seguir verificamos que a integral da função é 1 e fazemos o seu gráfico mostrado na Figura 18.
> 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
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.
> ef2 <- function(x){ x * f2(x) } > integrate(ef2, 0, 3)
1.333333 with absolute error < 7.3e-05
> 30 * integrate(ef2, 0, 3)$value
[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 19 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)
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.
(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} & \mbox { se $0 \leq x \leq 2$} \cr -\frac {3}{40}x + \frac {9}{20} & \mbox { se $ 2 < x \leq 6$} \cr 0 & \mbox { se $x < 0$ ou $x > 6$} \end {array} \right . \]
mostre que \(f(x)\) é uma f.d.p..
calcule os quartis da distribuição.
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.
qual a renda média nesta localidade?