Não foi possível enviar o arquivo. Será algum problema com as permissões?
Diferenças
Aqui você vê as diferenças entre duas revisões dessa página.
Ambos lados da revisão anterior Revisão anterior Próxima revisão | Revisão anterior Próxima revisão Ambos lados da revisão seguinte | ||
disciplinas:ce223:comandos2008 [2008/03/26 17:42] paulojus |
disciplinas:ce223:comandos2008 [2008/03/31 10:47] ehlers |
||
---|---|---|---|
Linha 545: | Linha 545: | ||
=== 24/03/2008 === | === 24/03/2008 === | ||
+ | **Problema:** | ||
+ | Considere duas variáveis aleatórias: ''N'' com distribuição Poisson de parâmetro 15 e ''Pr'' com distribuição exponencial de parâmetro 1/50000. Queremos obter o quantil 99 de uma variável aleatória ''Y'' dada pelo produto das anteriores, isto é ''Y = N * Pr''. | ||
+ | **Solução:** a solução analítica requer a obtenção da f.d.p. de ''Y''seguida da resolução de uma equação envolvendo integração desta f.d.p. para encontrar o quantil pedido. Entretanto, vamos aqui ilustrar uma alternativa numérica para resolver o problema, usando uma aproximação numérica por simulação. | ||
<code R> | <code R> | ||
N<- rpois(10000, lambda=15) | N<- rpois(10000, lambda=15) | ||
Linha 551: | Linha 554: | ||
Y<- N*Pr | Y<- N*Pr | ||
q99 <- quantile(Y, prob=0.99) | q99 <- quantile(Y, prob=0.99) | ||
- | hist(Y) | + | </code> |
+ | Vamos agora visualizar a distribuição de interesse de diferentes formas: pelo histograma das simulações e, | ||
+ | uma forma alternativa (e mais interessante!!!) utilizando estimação de densidades. | ||
+ | <code R> | ||
+ | hist(Y) | ||
hist(Y, prob=T) | hist(Y, prob=T) | ||
lines(density(Y)) | lines(density(Y)) | ||
plot(density(Y)) | plot(density(Y)) | ||
abline(v=q99) | abline(v=q99) | ||
+ | </code> | ||
- | + | Note que funções podem retornar resultados e/ou gráficos. A função ''hist()'' é um exemplo de função que retorna ambos. | |
- | hy<- hist(Y) | + | <code R> |
+ | hy <- hist(Y) | ||
hy | hy | ||
class(hy) | class(hy) | ||
Linha 567: | Linha 575: | ||
</code> | </code> | ||
- | Criando função: | + | Criando uma função -- um exemplo. Vamos encapsular todo o procedimento acima em uma função. Isto pode |
+ | ser útil para tornar a execução mais rápida e eficiente quando o procedimento deve ser repetido várias vezes. | ||
+ | (o equivalente a construir ''macros''). | ||
<code R> | <code R> | ||
mf<- function(n,lam,r,q){ | mf<- function(n,lam,r,q){ | ||
- | N<-rpois(n,lambda=lam) | + | N<-rpois(n,lambda=lam) |
- | Pr<-rexp(n,rate=r) | + | Pr<-rexp(n,rate=r) |
- | Y<-N*Pr | + | Y<-N*Pr |
- | qq <- quantile(Y,prob=q) | + | qq <- quantile(Y,prob=q) |
- | hist(Y,prob=T) | + | hist(Y,prob=T) |
- | lines(density(Y)) | + | lines(density(Y)) |
- | abline(v=qq) | + | abline(v=qq) |
- | text("topright", paste("quantil", q, "=", qq)) | + | text("topright", paste("quantil", q, "=", qq)) |
- | return(invisible(qq)) | + | return(invisible(qq)) |
} | } | ||
- | |||
mf(10000, 15, 1/50000, 0.99) | mf(10000, 15, 1/50000, 0.99) | ||
- | resp<- mf(10000,15, 1/50000, 0.99) | + | resp <- mf(10000,15, 1/50000, 0.99) |
resp | resp | ||
</code> | </code> | ||
+ | |||
+ | === 26/03/2008 === | ||
+ | Exercício proposto no material do curso e extensões discutidas em aula. | ||
+ | |||
+ | Calculando o valor da expressão | ||
+ | <code R> | ||
+ | x <- c(12, 11,14,15,10,11,14,11) | ||
+ | E=-length(x)*10 + sum(x) * log(10) - sum(log(factorial(x))) | ||
+ | E | ||
+ | </code> | ||
+ | |||
+ | Agora tornando mais flexível, escrevendo uma função que permite entrar com diferentes amostras e/ou valores de λ. | ||
+ | <code R> | ||
+ | mf <- function(y, lam){ | ||
+ | -length(y)*lam + sum(y) * log(lam) - sum(log(factorial(y))) | ||
+ | } | ||
+ | mf(y=x, lam=10) | ||
+ | </code> | ||
+ | |||
+ | Noque que está é a expressão da log-verossimilhanca para uma a.a. de uma distribuição de Poisson | ||
+ | <code R> | ||
+ | mf(y=x, lam=11) | ||
+ | mf(y=x, lam=12) | ||
+ | mf(y=x, lam=13) | ||
+ | </code> | ||
+ | |||
+ | Vamos então fazer o gráfico da função log-verossimilhança. Como esta é uma função do parâmetro λ vamos primeiro redefinir a função de uma forma mais conveniente colocando o parâmetro como o primeiro argumento. | ||
+ | <code R> | ||
+ | mf <- function(lam, y){ | ||
+ | -length(y)*lam + sum(y) * log(lam) - sum(log(factorial(y))) | ||
+ | } | ||
+ | l <- seq(5, 25, l=200) | ||
+ | ll <- mf(l) | ||
+ | ll <- mf(l, y=x) | ||
+ | plot(l, ll, type="l", xlab=expression(lambda), ylab=expression(l(lambda))) | ||
+ | </code> | ||
+ | |||
+ | Vamos agora indicar a solução analítica. | ||
+ | <code R> | ||
+ | mean(x) | ||
+ | abline(v=mean(x)) | ||
+ | </code> | ||
+ | |||
+ | A solução também poderia ser obtida por otimização numérica. Isto não é vantajoso para este problema mas pode ser a solução em casos onde a solução analítica não é disponível. | ||
+ | <code> | ||
+ | optimize(mf, c(min(x), max(x)), maximum=T, y=x) | ||
+ | </code> | ||
+ |