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/04 17:09] paulojus |
disciplinas:ce223:comandos2008 [2008/03/26 23:20] paulojus |
||
---|---|---|---|
Linha 111: | Linha 111: | ||
</code> | </code> | ||
- | ==== Semana 1 ==== | + | ==== Semana 2 ==== |
=== 03/03/2008 === | === 03/03/2008 === | ||
Linha 162: | Linha 162: | ||
</code> | </code> | ||
+ | Argumentos de funções: ordem dos argumentos, nomes dos argumentos, casamento parcial de nomes, argumentos com //default// | ||
<code R> | <code R> | ||
+ | args(dbinom) | ||
+ | dbinom(x, 10, 0.3) | ||
+ | dbinom(x=x, size=10, prob=0.3) | ||
+ | dbinom(prob=0.3, x=x, size=10) | ||
+ | dbinom(size=10, prob=0.3, x=x) | ||
+ | dbinom(x, s=10, p=0.3) | ||
+ | dbinom(x, p=0.3, s=10) | ||
+ | dbinom(x, p=0.3, s=10) | ||
+ | dbinom(x, 10, 0.3, log=F) | ||
+ | dbinom(x, 10, 0.3, log=T) | ||
</code> | </code> | ||
+ | === 05/03/2008 === | ||
+ | |||
+ | Verificando máximos, mínimos e suas posições em um vetor.<code R> | ||
+ | pesos <- c(67, 83, 56, 91, 58, 47, 82, 75) | ||
+ | max(pesos) | ||
+ | min(pesos) | ||
+ | which(pesos == max(pesos)) | ||
+ | which.max(pesos) | ||
+ | which(pesos == min(pesos)) | ||
+ | which.min(pesos) | ||
+ | </code> | ||
+ | |||
+ | Testando pela ocorrência de ''NA''. Note o uso do caracter de negação ''!'' <code R> | ||
+ | dat <- c(43, 56, NA, 23, 48, 33, NA, 29, 33, 39) | ||
+ | is.na(dat) | ||
+ | !is.na(dat) | ||
+ | which(is.na(dat)) | ||
+ | which(!is.na(dat)) | ||
+ | dat1 <- dat[!is.na(dat)] | ||
+ | </code> | ||
+ | |||
+ | ==== Semana 3 ==== | ||
+ | === 10/03/2008 === | ||
+ | |||
+ | Outras funções de probabilidades<code R> | ||
+ | args(dchisq) | ||
+ | args(dt) | ||
+ | </code> | ||
+ | |||
+ | Algumas contantes e operadores<code R> | ||
+ | exp(1) | ||
+ | exp(3) | ||
+ | print(pi, dig=12) | ||
+ | print(pi, dig=12) | ||
+ | pi | ||
+ | options(digits=12) | ||
+ | pi | ||
+ | exp(1) | ||
+ | </code> | ||
+ | |||
+ | Examinar opções de ''options()''<code R> | ||
+ | options() | ||
+ | </code> | ||
+ | |||
+ | |||
+ | === 12/03/2008 === | ||
+ | |||
+ | Funções para verificar exitência de ''NA'''s ''NaN'''s e valores infinitos (''Inf'') <code R> | ||
+ | x <- c(5, 0, -2) | ||
+ | x/0 | ||
+ | is.nan(x/0) | ||
+ | is.finite(x/0) | ||
+ | !is.finite(x/0) | ||
+ | </code> | ||
+ | |||
+ | Explicações sobre como o R armazana objetos (RAM e/ou dispositivos como por exemplo o HD)<code R> | ||
+ | save.image() | ||
+ | q() | ||
+ | </code> | ||
+ | |||
+ | Listando e apagando objetos | ||
<code R> | <code R> | ||
+ | ls() | ||
+ | objects() | ||
+ | rm(x) | ||
+ | rm(list=ls()) | ||
+ | args(ls) | ||
+ | ls() | ||
+ | ls(all=T) | ||
+ | </code> | ||
+ | |||
+ | Tipos de objetos no R: matrizes | ||
+ | <code R> | ||
+ | x <- 1:24 | ||
+ | x | ||
+ | m <- matrix(x, nr=6) | ||
+ | m | ||
+ | m <- matrix(x, nr=6, byrow=T) | ||
+ | m | ||
+ | x <- 1:25 | ||
+ | m <- matrix(x, nr=6) | ||
+ | m | ||
+ | attributes(x) | ||
+ | attributes(m) | ||
+ | </code> | ||
+ | |||
+ | Tipos de objetos no R: arrays | ||
+ | <code R> | ||
+ | x <- 1:24 | ||
+ | a <- array(x, dim=c(4,3,2)) | ||
+ | a | ||
+ | </code> | ||
+ | |||
+ | Digitação e conversão de uma tabela de tripla entrada (dada no quadro durante a aula) em um objeto do tipo ''array'' | ||
+ | | | PR || SC || RS || | ||
+ | | | Masculino | Feminino | Masculino | Feminino | Masculino | Feminino | | ||
+ | |Não Fuma ^ 45 ^ 16 ^ 21 ^ 33 ^ 40 ^ 45 ^ | ||
+ | |Fuma pouco ^ 28 ^ 22 ^ 34 ^ 21 ^ 50 ^ 37 ^ | ||
+ | |Fuma muito ^ 37 ^ 15 ^ 56 ^ 30 ^ 85 ^ 29 ^ | ||
+ | |||
+ | Comentários sobre ordem de entrada dos dados, cliclagem das variáveis e definição das dimensões do array<code R> | ||
+ | freqs <- scan() | ||
+ | 1: 45 | ||
+ | 2: 28 | ||
+ | 3: 37 | ||
+ | 4: 16 | ||
+ | 5: 22 | ||
+ | 6: 15 | ||
+ | 7: 21 | ||
+ | 8: 34 | ||
+ | 9: 56 | ||
+ | 10: 33 | ||
+ | 11: 21 | ||
+ | 12: 30 | ||
+ | 13: 40 | ||
+ | 14: 50 | ||
+ | 15: 85 | ||
+ | 16: 45 | ||
+ | 17: 37 | ||
+ | 18: 29 | ||
+ | 19: | ||
+ | freqs | ||
+ | Af <- array(freqs, dim=c(3,2,3)) | ||
+ | Af | ||
+ | </code> | ||
+ | |||
+ | ==== Semana 4 ==== | ||
+ | === 17/03/2008 === | ||
+ | <code R> | ||
+ | freqs = scan(file='http://leg.ufpr.br/~ehlers/CE223/fumo.dat') | ||
+ | |||
+ | array(freqs, dim=c(3,2,3)) | ||
+ | |||
+ | nomes = list(c('PR','SC','RS'), c('M','F'), c('nao fuma','fuma pouco','fuma muito')) | ||
+ | |||
+ | hf = array(freqs, dim=c(3,2,3), dimnames=nomes) | ||
+ | |||
+ | hf | ||
+ | |||
+ | m1 <- matrix(1:12, ncol = 3) | ||
+ | m1 | ||
+ | |||
+ | dimnames(m1) | ||
+ | |||
+ | dimnames(m1) <- list(c("L1", "L2", "L3", "L4"), c("C1", "C2", "C3")) | ||
+ | |||
+ | m1 | ||
+ | |||
+ | m2 <- cbind(1:5, 6:10) | ||
+ | m2 | ||
+ | |||
+ | m3 <- cbind(1:5, 6) | ||
+ | m3 | ||
+ | |||
+ | margin.table(m1, margin = 1) | ||
+ | |||
+ | apply(m1, 1, sum) | ||
+ | |||
+ | rowSums(m1) | ||
+ | |||
+ | margin.table(m1, margin = 2) | ||
+ | |||
+ | apply(m1, 2, sum) | ||
+ | |||
+ | colSums(m1) | ||
+ | </code> | ||
+ | |||
+ | Operacoes com matrizes | ||
+ | <code R> | ||
+ | m4 <- matrix(1:6, nc = 3) | ||
+ | m5 <- matrix(10 * (1:6), nc = 3) | ||
+ | m4 | ||
+ | |||
+ | m5 | ||
+ | |||
+ | m4 + m5 | ||
+ | |||
+ | m4 * m5 | ||
+ | |||
+ | m5 - m4 | ||
+ | |||
+ | m5/m4 | ||
+ | |||
+ | m4 %*% m5 | ||
+ | |||
+ | t(m4) | ||
+ | |||
+ | m6 = t(m4)%*% m5 | ||
+ | |||
+ | solve(m6) | ||
+ | |||
+ | m6[3,3]=20 | ||
+ | |||
+ | solve(m6) | ||
+ | |||
+ | mat <- matrix(c(1, 5, 2, 3, -2, 1, -1, 1, -1), nc = 3) | ||
+ | |||
+ | vec <- c(10, 15, 7) | ||
+ | |||
+ | solve(mat, vec) | ||
+ | </code> | ||
+ | |||
+ | === 19/03/2008 === | ||
+ | Data frames | ||
+ | <code R> | ||
+ | d1 = data.frame(x=1:10,y=c(51,54,61,67,68,75,77,75,80,82)) | ||
+ | |||
+ | d1 | ||
+ | |||
+ | x y | ||
+ | 1 1 51 | ||
+ | 2 2 54 | ||
+ | 3 3 61 | ||
+ | 4 4 67 | ||
+ | 5 5 68 | ||
+ | 6 6 75 | ||
+ | 7 7 77 | ||
+ | 8 8 75 | ||
+ | 9 9 80 | ||
+ | 10 10 82 | ||
+ | |||
+ | names(d1) | ||
+ | |||
+ | d1$x | ||
+ | d1$y | ||
+ | |||
+ | d1[,1] | ||
+ | d1[,2] | ||
+ | |||
+ | plot(d1) | ||
+ | |||
+ | plot(d1$x,d1$y) | ||
+ | |||
+ | d2 = data.frame(Y = c(rnorm(5,mean=10,sd=2), rnorm(5,16,2), rnorm(5,14,2))) | ||
+ | |||
+ | gl(3, 5) | ||
+ | [1] 1 1 1 1 1 2 2 2 2 2 3 3 3 3 3 | ||
+ | Levels: 1 2 3 | ||
+ | |||
+ | d2$lev = gl(3, 5) | ||
+ | |||
+ | d2 | ||
+ | |||
+ | is.factor(d2$lev) | ||
+ | </code> | ||
+ | |||
+ | Aplicando uma funcao a um Data Frame separando por fatores | ||
+ | <code R> | ||
+ | by(d2$Y, d2$lev, summary) | ||
+ | </code> | ||
+ | |||
+ | Criando um Data Frame a partir de todas as combinacoes de 2 fatores. | ||
+ | <code R> | ||
+ | d3 = expand.grid(1:3, 4:5) | ||
+ | d3 | ||
+ | |||
+ | is.data.frame(d3) | ||
+ | </code> | ||
+ | |||
+ | Adicionando colunas | ||
+ | <code R> | ||
+ | d4 = data.frame(peso=rnorm(15,65,5),altura=rnorm(15,160,10)) | ||
+ | d4 | ||
+ | |||
+ | d4=cbind(d4,sexo=c(rep('M',10),rep('F',5))) | ||
+ | d4 | ||
+ | </code> | ||
+ | |||
+ | Toda coluna que não seja composta exclusivamente de números é definida como um fator. | ||
+ | <code R> | ||
+ | is.factor(d4$sexo) | ||
+ | [1] TRUE | ||
+ | |||
+ | letters | ||
+ | |||
+ | LETTERS | ||
+ | |||
+ | d4=cbind(d4,nome=letters[1:15]) | ||
+ | |||
+ | is.factor(d4$nome) | ||
+ | [1] TRUE | ||
+ | |||
+ | d4$nome= as.character(d4$nome) | ||
+ | |||
+ | d4 | ||
+ | |||
+ | dim(d4) | ||
+ | |||
+ | names(d4) | ||
+ | |||
+ | dimnames(d4) | ||
+ | |||
+ | rownames(d4) | ||
+ | |||
+ | colnames(d4) | ||
+ | |||
+ | d4[d4$sexo=='M',1:2] | ||
+ | |||
+ | d4[d4$sexo=='F',4] | ||
+ | |||
+ | by(d4[,1:2],d4$sexo,function(x)x) | ||
+ | |||
+ | by(d4[,4],d4$sexo,function(x)as.character(x)) | ||
+ | </code> | ||
+ | |||
+ | Listas | ||
+ | <code R> | ||
+ | |||
+ | Listas sao estruturas genericas e flexiveis que permitem armazenar | ||
+ | diversos formatos em um unico objeto. | ||
+ | |||
+ | lis1 <- list(A = 1:10, B = "CE 223", C = matrix(1:9,ncol = 3)) | ||
+ | lis1 | ||
+ | |||
+ | $A | ||
+ | [1] 1 2 3 4 5 6 7 8 9 10 | ||
+ | |||
+ | $B | ||
+ | [1] "CE 223" | ||
+ | | ||
+ | $C | ||
+ | [,1] [,2] [,3] | ||
+ | [1,] 1 4 7 | ||
+ | [2,] 2 5 8 | ||
+ | [3,] 3 6 9 | ||
+ | </code> | ||
+ | |||
+ | Varias funcoes do R retornam listas | ||
+ | <code R> | ||
+ | d1 | ||
+ | |||
+ | lis2 = lm (y ~ x, data=d1) | ||
+ | |||
+ | lis2 | ||
+ | |||
+ | is.list(lis2) | ||
+ | |||
+ | class(lis2) | ||
+ | |||
+ | summary(lis2) | ||
+ | |||
+ | anova(lis2) | ||
+ | |||
+ | names(lis2) | ||
+ | |||
+ | lis2$pred | ||
+ | |||
+ | lis2$residuals | ||
+ | |||
+ | par(mfrow=c(2,2)) | ||
+ | |||
+ | plot(lis2) | ||
+ | </code> | ||
+ | |||
+ | Selecionando elementos de uma lista | ||
+ | <code R> | ||
+ | lis1$A | ||
+ | |||
+ | lis2$coeff | ||
+ | |||
+ | lis1[3] | ||
+ | |||
+ | lis1[[3]] | ||
+ | |||
+ | </code> | ||
+ | |||
+ | ==== Semana 5 ==== | ||
+ | |||
+ | === 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> | ||
+ | N<- rpois(10000, lambda=15) | ||
+ | Pr<- rexp(10000, rate=1/50000) | ||
+ | Y<- N*Pr | ||
+ | q99 <- quantile(Y, prob=0.99) | ||
+ | </code> | ||
+ | |||
+ | Vamosagora visualizar a distribuição de interesse de diferentes formas: pelo hiostograma das simulações e, | ||
+ | uma forma alternativa (e mais interessante!!!) utilizando estimação de densidades. | ||
+ | <code R> | ||
+ | hist(Y) | ||
+ | hist(Y, prob=T) | ||
+ | lines(density(Y)) | ||
+ | plot(density(Y)) | ||
+ | abline(v=q99) | ||
+ | </code> | ||
+ | |||
+ | Note que que funçãos podem retornar resultados e/ou gráficos. A função ''hist()'' é um exemplo de função que retorna ambos. | ||
+ | <code R> | ||
+ | hy<- hist(Y) | ||
+ | hy | ||
+ | class(hy) | ||
+ | plot(hy) | ||
+ | hist | ||
+ | </code> | ||
+ | |||
+ | Criando uma função-- um exemplo. Vamosencapsular todo o procedimento acima em uma função. Isto pode | ||
+ | ser útil para tornar a execução mias rápida e eficiente quando o procedimento deve ser repetido várias vezes. | ||
+ | (o equivalente a construir ''macros''). | ||
+ | <code R> | ||
+ | mf<- function(n,lam,r,q){ | ||
+ | N<-rpois(n,lambda=lam) | ||
+ | Pr<-rexp(n,rate=r) | ||
+ | Y<-N*Pr | ||
+ | qq <- quantile(Y,prob=q) | ||
+ | hist(Y,prob=T) | ||
+ | lines(density(Y)) | ||
+ | abline(v=qq) | ||
+ | text("topright", paste("quantil", q, "=", qq)) | ||
+ | return(invisible(qq)) | ||
+ | } | ||
+ | mf(10000, 15, 1/50000, 0.99) | ||
+ | resp <- mf(10000,15, 1/50000, 0.99) | ||
+ | resp | ||
+ | </code> | ||
+ | |||
+ | === 24/03/2008 === | ||
+ | Exercício proposto no material do cursoe 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 cem casosonde asolução analítica não é disponível. | ||
+ | <code> | ||
+ | optimize(mf, c(min(x), max(x)), maximum=T, y=x) | ||
</code> | </code> | ||