Curso de estatística experimental com aplicações em R |
|
|
12 à 14 de Novembro de 2014 - Manaus - AM |
Prof. Dr. Walmes M. Zeviani |
|
Embrapa Amazônia Ocidental |
Lab. de Estatística e Geoinformação - LEG |
| Departamento de Estatística - UFPR | |
##-----------------------------------------------------------------------------
## Sequências regulares.
1:5## [1] 1 2 3 4 52*(1:5)-1## [1] 1 3 5 7 93*(1:5)-2## [1] 1 4 7 10 13x <- seq(from=1, to=12, by=2); x## [1] 1 3 5 7 9 11x <- seq(from=1, to=12, length.out=5); x## [1] 1.00 3.75 6.50 9.25 12.00x <- seq(from=1, by=3, length.out=5); x## [1] 1 4 7 10 13x <- seq(to=20, by=3, length.out=6); x## [1] 5 8 11 14 17 20##-----------------------------------------------------------------------------
## Repetições regulares.
x <- rep(1, times=5); x## [1] 1 1 1 1 1x <- rep(1:3, times=3); x## [1] 1 2 3 1 2 3 1 2 3x <- rep(1:3, each=3); x## [1] 1 1 1 2 2 2 3 3 3##-----------------------------------------------------------------------------
## Grides.
## da <- expand.grid(x1=1:3, x2=1:3)
da <- expand.grid(##bloco=1:3,
vari=c("Carmen","Lucia","Teresa"),
## bloco=factor(1:3),
bloco=gl(3,1),
nitro=seq(0, 100, by=20))
da## vari bloco nitro
## 1 Carmen 1 0
## 2 Lucia 1 0
## 3 Teresa 1 0
## 4 Carmen 2 0
## 5 Lucia 2 0
## 6 Teresa 2 0
## 7 Carmen 3 0
## 8 Lucia 3 0
## 9 Teresa 3 0
## 10 Carmen 1 20
## 11 Lucia 1 20
## 12 Teresa 1 20
## 13 Carmen 2 20
## 14 Lucia 2 20
## 15 Teresa 2 20
## 16 Carmen 3 20
## 17 Lucia 3 20
## 18 Teresa 3 20
## 19 Carmen 1 40
## 20 Lucia 1 40
## 21 Teresa 1 40
## 22 Carmen 2 40
## 23 Lucia 2 40
## 24 Teresa 2 40
## 25 Carmen 3 40
## 26 Lucia 3 40
## 27 Teresa 3 40
## 28 Carmen 1 60
## 29 Lucia 1 60
## 30 Teresa 1 60
## 31 Carmen 2 60
## 32 Lucia 2 60
## 33 Teresa 2 60
## 34 Carmen 3 60
## 35 Lucia 3 60
## 36 Teresa 3 60
## 37 Carmen 1 80
## 38 Lucia 1 80
## 39 Teresa 1 80
## 40 Carmen 2 80
## 41 Lucia 2 80
## 42 Teresa 2 80
## 43 Carmen 3 80
## 44 Lucia 3 80
## 45 Teresa 3 80
## 46 Carmen 1 100
## 47 Lucia 1 100
## 48 Teresa 1 100
## 49 Carmen 2 100
## 50 Lucia 2 100
## 51 Teresa 2 100
## 52 Carmen 3 100
## 53 Lucia 3 100
## 54 Teresa 3 100str(da)## 'data.frame': 54 obs. of 3 variables:
## $ vari : Factor w/ 3 levels "Carmen","Lucia",..: 1 2 3 1 2 3 1 2 3 1 ...
## $ bloco: Factor w/ 3 levels "1","2","3": 1 1 1 2 2 2 3 3 3 1 ...
## $ nitro: num 0 0 0 0 0 0 0 0 0 20 ...
## - attr(*, "out.attrs")=List of 2
## ..$ dim : Named int 3 3 6
## .. ..- attr(*, "names")= chr "vari" "bloco" "nitro"
## ..$ dimnames:List of 3
## .. ..$ vari : chr "vari=Carmen" "vari=Lucia" "vari=Teresa"
## .. ..$ bloco: chr "bloco=1" "bloco=2" "bloco=3"
## .. ..$ nitro: chr "nitro= 0" "nitro= 20" "nitro= 40" "nitro= 60" ...da <- expand.grid(trat=c("tratado","controle"),
bloco=c("I","II","III"))
da## trat bloco
## 1 tratado I
## 2 controle I
## 3 tratado II
## 4 controle II
## 5 tratado III
## 6 controle III##-----------------------------------------------------------------------------
## Números uniformes.
x <- runif(10); x## [1] 0.005343402 0.163498536 0.468872580 0.896007781 0.026639577 0.550113728 0.219849132
## [8] 0.300428519 0.901197906 0.656097755##-----------------------------------------------------------------------------
## Números da distribuição normal padrão.
x <- rnorm(10, mean=0, sd=1)
## curve(dnorm(x, 0, 1), -4, 4)
## rug(x)
##-----------------------------------------------------------------------------
## Amostra aleatória de população discreta finita com e sem reposição.
u <- c("Superman","Batman","Aquaman","Flash","Green Lantern")
sample(u, size=3)## [1] "Batman" "Green Lantern" "Aquaman"u <- c("Iron Man","Thor","Captain America","Hulk","Black Widow","Hawkeye")
sample(u, size=10, replace=TRUE)## [1] "Black Widow" "Iron Man" "Hulk" "Captain America"
## [5] "Black Widow" "Hulk" "Captain America" "Iron Man"
## [9] "Black Widow" "Hawkeye"##-----------------------------------------------------------------------------
## Um vetor com número uniformes.
x <- runif(6)
## As posições do menor para o maior valor.
order(x)## [1] 1 4 5 6 3 2## O vetor ordenado do menor para o maior.
sort(x, decreasing=FALSE)## [1] 0.4622079 0.4827947 0.7711735 0.9045427 0.9196033 0.9350125x[order(x)]## [1] 0.4622079 0.4827947 0.7711735 0.9045427 0.9196033 0.9350125## Os valores em ordem inversa.
rev(x)## [1] 0.9045427 0.7711735 0.4827947 0.9196033 0.9350125 0.4622079## Ordenando um vetor por valores em outro.
u <- c("Iron Man","Thor","Captain America","Hulk","Black Widow","Hawkeye")
u[order(x)]## [1] "Iron Man" "Hulk" "Black Widow" "Hawkeye"
## [5] "Captain America" "Thor"## O vetor com os elementos aleatorizados.
sample(u)## [1] "Black Widow" "Hawkeye" "Thor" "Captain America"
## [5] "Iron Man" "Hulk"##-----------------------------------------------------------------------------
## Um vetor com dados de precipitação.
precip## Mobile Juneau Phoenix Little Rock
## 67.0 54.7 7.0 48.5
## Los Angeles Sacramento San Francisco Denver
## 14.0 17.2 20.7 13.0
## Hartford Wilmington Washington Jacksonville
## 43.4 40.2 38.9 54.5
## Miami Atlanta Honolulu Boise
## 59.8 48.3 22.9 11.5
## Chicago Peoria Indianapolis Des Moines
## 34.4 35.1 38.7 30.8
## Wichita Louisville New Orleans Portland
## 30.6 43.1 56.8 40.8
## Baltimore Boston Detroit Sault Ste. Marie
## 41.8 42.5 31.0 31.7
## Duluth Minneapolis/St Paul Jackson Kansas City
## 30.2 25.9 49.2 37.0
## St Louis Great Falls Omaha Reno
## 35.9 15.0 30.2 7.2
## Concord Atlantic City Albuquerque Albany
## 36.2 45.5 7.8 33.4
## Buffalo New York Charlotte Raleigh
## 36.1 40.2 42.7 42.5
## Bismark Cincinati Cleveland Columbus
## 16.2 39.0 35.0 37.0
## Oklahoma City Portland Philadelphia Pittsburg
## 31.4 37.6 39.9 36.2
## Providence Columbia Sioux Falls Memphis
## 42.8 46.4 24.7 49.1
## Nashville Dallas El Paso Houston
## 46.0 35.9 7.8 48.2
## Salt Lake City Burlington Norfolk Richmond
## 15.2 32.5 44.7 42.6
## Seattle Tacoma Spokane Charleston Milwaukee
## 38.8 17.4 40.8 29.1
## Cheyenne San Juan
## 14.6 59.2x <- precip
sum(x)## [1] 2442mean(x)## [1] 34.88571median(x)## [1] 36.6min(x)## [1] 7max(x)## [1] 67range(x)## [1] 7 67fivenum(x)## Phoenix Milwaukee Pittsburg Providence Mobile
## 7.0 29.1 36.6 42.8 67.0IQR(x)## [1] 13.4quantile(x, probs=c(0.05, 0.95))## 5% 95%
## 9.465 55.855var(x)## [1] 187.8723sd(x)## [1] 13.70665mad(x)## [1] 9.56277diff(range(x))## [1] 60100*sd(x)/mean(x)## [1] 39.29015require(fBasics)## Loading required package: fBasics
## Loading required package: timeDate
## Loading required package: methods
## Loading required package: timeSeries
##
##
## Rmetrics Package fBasics
## Analysing Markets and calculating Basic Statistics
## Copyright (C) 2005-2014 Rmetrics Association Zurich
## Educational Software for Financial Engineering and Computational Science
## Rmetrics is free software and comes with ABSOLUTELY NO WARRANTY.
## https://www.rmetrics.org --- Mail to: info@rmetrics.orgskewness(x)## [1] -0.2852747
## attr(,"method")
## [1] "moment"kurtosis(x)## [1] -0.38499
## attr(,"method")
## [1] "excess"require(EnvStats)## Loading required package: EnvStats
##
## Attaching package: 'EnvStats'
##
## The following objects are masked from 'package:timeDate':
##
## kurtosis, skewness
##
## The following object is masked from 'package:stats':
##
## predict.lm## Algumas funções são substituídas.
skewness(x)## [1] -0.2979212kurtosis(x)## [1] -0.2410105cv(x)## [1] 0.3929015summaryStats(x)## N Mean SD Median Min Max
## x 70 34.8857 13.7067 36.6 7 67sx <- summaryFull(x)
sx## x
## N 70
## Mean 34.89
## Median 36.6
## 10% Trimmed Mean 35.22
## Geometric Mean 31.26
## Skew -0.2979
## Kurtosis -0.241
## Min 7
## Max 67
## Range 60
## 1st Quartile 29.38
## 3rd Quartile 42.77
## Standard Deviation 13.71
## Geometric Standard Deviation 1.696
## Interquartile Range 13.39
## Median Absolute Deviation 9.563
## Coefficient of Variation 0.3929str(sx)## summaryStats [1:17, 1] 70 34.9 36.6 35.2 31.3 ...
## - attr(*, "dimnames")=List of 2
## ..$ : chr [1:17] "N" "Mean" "Median" "10% Trimmed Mean" ...
## ..$ : chr "x"
## - attr(*, "stats.in.rows")= logi TRUE
## - attr(*, "drop0trailing")= logi TRUEsx["10% Trimmed Mean",]## [1] 35.22A família *apply e agregados representa um conjunto de funções destinadas à tarefas por estrato e/ou margem (índice). Tarefas assim foram recentemente batizadas de SAC (split-apply-combine). São essas funções do pacote básico do R que permiter obter médias por tratamento, coeficientes de variação por resposta (por coluna). Isoladamente cada uma delas têm seu papel e juntas resolvem uma série de problemas.
##-----------------------------------------------------------------------------
## Os membros da família *apply.
## De uso simples e frequente.
## apply : on Arrays margins (não é aaply para evitar cacofonia).
## lapply : on Lists itens or vector elements or data.frame columns.
## sapply : lapply that Simplifies when possible.
## tapply : on a ragged array, return as Tabular format.
## De uso menos frequente.
## mapply : lapply over Multiple list or vector arguments.
## eapply : on Environments.
## rapply : Recursive.
## vapply : Vectorized.
## dendrapply : related to dendrogramns.
## kernapply : related kernel.
## Além destas, tem-se outras funções que são os agregados da família.
## by(): divide de acordo com estrato e aplica algo.
## aggragate(): divide por estrato e aplica por coluna.
## sweep(): operações entre matrizes e vetores.
##-----------------------------------------------------------------------------
## tapply.
str(npk)## 'data.frame': 24 obs. of 5 variables:
## $ block: Factor w/ 6 levels "1","2","3","4",..: 1 1 1 1 2 2 2 2 3 3 ...
## $ N : Factor w/ 2 levels "0","1": 1 2 1 2 2 2 1 1 1 2 ...
## $ P : Factor w/ 2 levels "0","1": 2 2 1 1 1 2 1 2 2 2 ...
## $ K : Factor w/ 2 levels "0","1": 2 1 1 2 1 2 2 1 1 2 ...
## $ yield: num 49.5 62.8 46.8 57 59.8 58.5 55.5 56 62.8 55.8 ...## Média de yield para cada nível de N.
r <- tapply(npk$yield, npk$N, mean, na.rm=TRUE)
r## 0 1
## 52.06667 57.68333class(r)## [1] "array"## Média de yield para combinando os níveis de N e P.
r <- tapply(npk$yield, list(npk$N, npk$P), mean)
r## 0 1
## 0 51.71667 52.41667
## 1 59.21667 56.15000## Melhor usa com with() para simplificar a declaração.
r <- with(npk, tapply(yield, list(N, P), mean))
r## 0 1
## 0 51.71667 52.41667
## 1 59.21667 56.15000## O mesmo considerando N, P e K.
r <- with(npk, tapply(yield, list(N, P, K), mean))
r## , , 0
##
## 0 1
## 0 51.43333 54.33333
## 1 63.76667 57.93333
##
## , , 1
##
## 0 1
## 0 52.00000 50.50000
## 1 54.66667 54.36667## Nomes na lista geram nomes para as dimensões do array.
r <- with(npk, tapply(yield, list(Nitro=N, Phos=P, Pot=K), mean))
r## , , Pot = 0
##
## Phos
## Nitro 0 1
## 0 51.43333 54.33333
## 1 63.76667 57.93333
##
## , , Pot = 1
##
## Phos
## Nitro 0 1
## 0 52.00000 50.50000
## 1 54.66667 54.36667##-----------------------------------------------------------------------------
## aggregate.
## A aggragate funciona com uso de formula, além de poder ser usada como a
## tapply(). O resultado é em data.frame.
s <- with(npk,
aggregate(cbind(Y=yield), list(Nitro=N, Phos=P, Pot=K), mean))
s## Nitro Phos Pot Y
## 1 0 0 0 51.43333
## 2 1 0 0 63.76667
## 3 0 1 0 54.33333
## 4 1 1 0 57.93333
## 5 0 0 1 52.00000
## 6 1 0 1 54.66667
## 7 0 1 1 50.50000
## 8 1 1 1 54.36667## Será obtido o mesmo, mas usando uma fórmula para representar o que se
## deseja.
s <- aggregate(yield~N+P+K, data=npk, mean)
s## N P K yield
## 1 0 0 0 51.43333
## 2 1 0 0 63.76667
## 3 0 1 0 54.33333
## 4 1 1 0 57.93333
## 5 0 0 1 52.00000
## 6 1 0 1 54.66667
## 7 0 1 1 50.50000
## 8 1 1 1 54.36667## Diferente da tapply, a aggregate pode ter mais de uma variáveis
## resposta. Por falta de outra variável resposta, será usando o log de
## yield.
s <- aggregate(cbind(y=yield, log.y=log(yield))~N+P+K, data=npk, mean)
s## N P K y log.y
## 1 0 0 0 51.43333 3.937606
## 2 1 0 0 63.76667 4.153156
## 3 0 1 0 54.33333 3.984677
## 4 1 1 0 57.93333 4.056245
## 5 0 0 1 52.00000 3.947143
## 6 1 0 1 54.66667 3.999207
## 7 0 1 1 50.50000 3.921254
## 8 1 1 1 54.36667 3.992844##-----------------------------------------------------------------------------
## by.
by(data=npk, INDICES=with(npk, N), FUN=nrow)## with(npk, N): 0
## [1] 12
## -------------------------------------------------------------------
## with(npk, N): 1
## [1] 12r <- with(npk, by(yield, INDICES=N, FUN=mean)); r## N: 0
## [1] 52.06667
## -------------------------------------------------------------------
## N: 1
## [1] 57.68333str(r)## by [1:2(1d)] 52.1 57.7
## - attr(*, "dimnames")=List of 1
## ..$ N: chr [1:2] "0" "1"
## - attr(*, "call")= language by.default(data = yield, INDICES = N, FUN = mean)c(is.list(r), is.array(r))## [1] FALSE TRUEr <- with(npk, by(yield, INDICES=list(N=N, P=P, K=K), FUN=mean)); r## N: 0
## P: 0
## K: 0
## [1] 51.43333
## -------------------------------------------------------------------
## N: 1
## P: 0
## K: 0
## [1] 63.76667
## -------------------------------------------------------------------
## N: 0
## P: 1
## K: 0
## [1] 54.33333
## -------------------------------------------------------------------
## N: 1
## P: 1
## K: 0
## [1] 57.93333
## -------------------------------------------------------------------
## N: 0
## P: 0
## K: 1
## [1] 52
## -------------------------------------------------------------------
## N: 1
## P: 0
## K: 1
## [1] 54.66667
## -------------------------------------------------------------------
## N: 0
## P: 1
## K: 1
## [1] 50.5
## -------------------------------------------------------------------
## N: 1
## P: 1
## K: 1
## [1] 54.36667str(r)## by [1:2, 1:2, 1:2] 51.4 63.8 54.3 57.9 52 ...
## - attr(*, "dimnames")=List of 3
## ..$ N: chr [1:2] "0" "1"
## ..$ P: chr [1:2] "0" "1"
## ..$ K: chr [1:2] "0" "1"
## - attr(*, "call")= language by.default(data = yield, INDICES = list(N = N, P = P, K = K), FUN = mean)c(is.list(r), is.array(r))## [1] FALSE TRUEclass(r) ## Como é de classe by ele é mostrado de forma diferente.## [1] "by"unclass(r) ## Se a classe é removida, então é mostrado como array comum.## , , K = 0
##
## P
## N 0 1
## 0 51.43333 54.33333
## 1 63.76667 57.93333
##
## , , K = 1
##
## P
## N 0 1
## 0 52.00000 50.50000
## 1 54.66667 54.36667
##
## attr(,"call")
## by.default(data = yield, INDICES = list(N = N, P = P, K = K),
## FUN = mean)Em resumo, as funções tapply, aggregate e by fazem tarefas por estrato. Em outras palavras, separam os valores em um (ou mais) vetores respeitando valores em outro (ou mais) e em seguida aplicam uma função. A diferença é como declarar e o que é retornado.
##-----------------------------------------------------------------------------
## apply.
Titanic## , , Age = Child, Survived = No
##
## Sex
## Class Male Female
## 1st 0 0
## 2nd 0 0
## 3rd 35 17
## Crew 0 0
##
## , , Age = Adult, Survived = No
##
## Sex
## Class Male Female
## 1st 118 4
## 2nd 154 13
## 3rd 387 89
## Crew 670 3
##
## , , Age = Child, Survived = Yes
##
## Sex
## Class Male Female
## 1st 5 1
## 2nd 11 13
## 3rd 13 14
## Crew 0 0
##
## , , Age = Adult, Survived = Yes
##
## Sex
## Class Male Female
## 1st 57 140
## 2nd 14 80
## 3rd 75 76
## Crew 192 20str(Titanic)## table [1:4, 1:2, 1:2, 1:2] 0 0 35 0 0 0 17 0 118 154 ...
## - attr(*, "dimnames")=List of 4
## ..$ Class : chr [1:4] "1st" "2nd" "3rd" "Crew"
## ..$ Sex : chr [1:2] "Male" "Female"
## ..$ Age : chr [1:2] "Child" "Adult"
## ..$ Survived: chr [1:2] "No" "Yes"is.array(Titanic)## [1] TRUEdimnames(Titanic)## $Class
## [1] "1st" "2nd" "3rd" "Crew"
##
## $Sex
## [1] "Male" "Female"
##
## $Age
## [1] "Child" "Adult"
##
## $Survived
## [1] "No" "Yes"sum(Titanic[ ,1, , ]) ## Total de homens.## [1] 1731sum(Titanic[ ,2, , ]) ## Total de mulheres.## [1] 470apply(Titanic, MARGIN=2, sum) ## Totais das margens para Sex.## Male Female
## 1731 470apply(Titanic, MARGIN=1, sum) ## Para Class.## 1st 2nd 3rd Crew
## 325 285 706 885apply(Titanic, MARGIN=c(1,2), sum) ## Class e Sex.## Sex
## Class Male Female
## 1st 180 145
## 2nd 179 106
## 3rd 510 196
## Crew 862 23apply(Titanic, MARGIN=c(2,4), sum) ## Sex e Survived.## Survived
## Sex No Yes
## Male 1364 367
## Female 126 344apply(Titanic, MARGIN=c(3,4), sum) ## Age e Survived.## Survived
## Age No Yes
## Child 52 57
## Adult 1438 654str(HairEyeColor)## table [1:4, 1:4, 1:2] 32 53 10 3 11 50 10 30 10 25 ...
## - attr(*, "dimnames")=List of 3
## ..$ Hair: chr [1:4] "Black" "Brown" "Red" "Blond"
## ..$ Eye : chr [1:4] "Brown" "Blue" "Hazel" "Green"
## ..$ Sex : chr [1:2] "Male" "Female"dimnames(HairEyeColor)## $Hair
## [1] "Black" "Brown" "Red" "Blond"
##
## $Eye
## [1] "Brown" "Blue" "Hazel" "Green"
##
## $Sex
## [1] "Male" "Female"apply(HairEyeColor, MARGIN=1, sum) ## Por cor de cabelo.## Black Brown Red Blond
## 108 286 71 127apply(HairEyeColor, MARGIN=2, sum) ## Por cor de olhos.## Brown Blue Hazel Green
## 220 215 93 64apply(HairEyeColor, MARGIN=3, sum) ## Por cor de sexo.## Male Female
## 279 313apply(HairEyeColor, MARGIN=c(1,2), sum) ## Por cor de cabelo.## Eye
## Hair Brown Blue Hazel Green
## Black 68 20 15 5
## Brown 119 84 54 29
## Red 26 17 14 14
## Blond 7 94 10 16##-----------------------------------------------------------------------------
## lapply e sapply.
is.list(rock)## [1] TRUEstr(rock) ## Todas as colunas tem conteúdo numérico.## 'data.frame': 48 obs. of 4 variables:
## $ area : int 4990 7002 7558 7352 7943 7979 9333 8209 8393 6425 ...
## $ peri : num 2792 3893 3931 3869 3949 ...
## $ shape: num 0.0903 0.1486 0.1833 0.1171 0.1224 ...
## $ perm : num 6.3 6.3 6.3 6.3 17.1 17.1 17.1 17.1 119 119 ...lapply(rock, mean) ## Média## $area
## [1] 7187.729
##
## $peri
## [1] 2682.212
##
## $shape
## [1] 0.2181104
##
## $perm
## [1] 415.45lapply(rock, range) ## Extremos.## $area
## [1] 1016 12212
##
## $peri
## [1] 308.642 4864.220
##
## $shape
## [1] 0.0903296 0.4641250
##
## $perm
## [1] 6.3 1300.0## Porque trata-se de um data.frame, dá pra usar apply também.
apply(rock, MARGIN=2, mean)## area peri shape perm
## 7187.7291667 2682.2119375 0.2181104 415.4500000apply(rock, MARGIN=2, range)## area peri shape perm
## [1,] 1016 308.642 0.0903296 6.3
## [2,] 12212 4864.220 0.4641250 1300.0sapply(rock, mean) ## Foi possível simplificar para um vetor.## area peri shape perm
## 7187.7291667 2682.2119375 0.2181104 415.4500000sapply(rock, range) ## Foi possível simplificar para uma matriz.## area peri shape perm
## [1,] 1016 308.642 0.0903296 6.3
## [2,] 12212 4864.220 0.4641250 1300.0str(iris)## 'data.frame': 150 obs. of 5 variables:
## $ Sepal.Length: num 5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 ...
## $ Sepal.Width : num 3.5 3 3.2 3.1 3.6 3.9 3.4 3.4 2.9 3.1 ...
## $ Petal.Length: num 1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 ...
## $ Petal.Width : num 0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 ...
## $ Species : Factor w/ 3 levels "setosa","versicolor",..: 1 1 1 1 1 1 1 1 1 1 ...sapply(iris, mean)## Warning in mean.default(X[[5L]], ...): argument is not numeric or logical: returning NA## Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 5.843333 3.057333 3.758000 1.199333 NAsapply(iris, summary)## $Sepal.Length
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 4.300 5.100 5.800 5.843 6.400 7.900
##
## $Sepal.Width
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 2.000 2.800 3.000 3.057 3.300 4.400
##
## $Petal.Length
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.000 1.600 4.350 3.758 5.100 6.900
##
## $Petal.Width
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.100 0.300 1.300 1.199 1.800 2.500
##
## $Species
## setosa versicolor virginica
## 50 50 50lapply(iris, is.numeric) ## Quais colunas tem conteúdo numérico?## $Sepal.Length
## [1] TRUE
##
## $Sepal.Width
## [1] TRUE
##
## $Petal.Length
## [1] TRUE
##
## $Petal.Width
## [1] TRUE
##
## $Species
## [1] FALSElapply(iris, class) ## Qual a classe?## $Sepal.Length
## [1] "numeric"
##
## $Sepal.Width
## [1] "numeric"
##
## $Petal.Length
## [1] "numeric"
##
## $Petal.Width
## [1] "numeric"
##
## $Species
## [1] "factor"## Também se pode usar apply pois iris é um data.frame.
apply(iris, MARGIN=2, is.numeric)## Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## FALSE FALSE FALSE FALSE FALSEsapply(iris, class) ## Foi possível simplificar.## Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## "numeric" "numeric" "numeric" "numeric" "factor"## Não foi possível simplificar.
sapply(iris, summary)## $Sepal.Length
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 4.300 5.100 5.800 5.843 6.400 7.900
##
## $Sepal.Width
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 2.000 2.800 3.000 3.057 3.300 4.400
##
## $Petal.Length
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.000 1.600 4.350 3.758 5.100 6.900
##
## $Petal.Width
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.100 0.300 1.300 1.199 1.800 2.500
##
## $Species
## setosa versicolor virginica
## 50 50 50## Separar as colunas que são númericas e então pedir o summary.
i <- sapply(iris, is.numeric); i## Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## TRUE TRUE TRUE TRUE FALSEsapply(iris[,i], summary) ## Simplificou.## Sepal.Length Sepal.Width Petal.Length Petal.Width
## Min. 4.300 2.000 1.000 0.100
## 1st Qu. 5.100 2.800 1.600 0.300
## Median 5.800 3.000 4.350 1.300
## Mean 5.843 3.057 3.758 1.199
## 3rd Qu. 6.400 3.300 5.100 1.800
## Max. 7.900 4.400 6.900 2.500Em resumo, as funções apply, lapply e sapply aplicam funções sobre índices, ou seja, em arranjos as operações são varrendo linhas, colunas, etc e em listas são varrendo os ítens. A diferença de lapply e sapply é que a última tenta simplificar o resultado acomodanto em um tipo de objeto mais simples que uma lista quando possível.
##-----------------------------------------------------------------------------
## Função para calcular a hipotenusa de um triângulo retângulo.
hipo <- function(a, b){ ## Argumentos.
h <- sqrt(a^2+b^2) ## Processo/procedimento.
return(h) ## Retorno.
}
## Explorando o objeto.
class(hipo)
mode(hipo)
args(hipo)
body(hipo)
str(hipo)
## Usando a função.
hipo(3, 4)
hipo(3, 4.6)
hipo(3, 3)
## Note que ela opera vetorialmente.
hipo(a=3, b=c(4, 4.6, 3))
##-----------------------------------------------------------------------------
## Abrindo algumas funções.
body(plot) ## Indica que é uma função genérica.
methods(generic.function="plot") ## As várias 'faces' da plot.
methods(generic.function="summary") ## As várias 'faces' da plot.
## Abrindo a função.
plot.default
plot.ecdf
## plot.ts
## As com asteriscos podem ser abertas assim.
getAnywhere(plot.histogram)
## getS3method(f="plot", class="histogram")
## Abrindo a função sequência.
methods(generic.function="seq")
getAnywhere(seq.default)
## Abrindo a matrix.
matrix
runif ## Chama uma função escrita em C.
##-----------------------------------------------------------------------------
## Incluindo o argumento ângulo.
hipo2 <- function(a, b, angle){
h2 <- a^2+b^2-2*cos(angle)
return(sqrt(h2))
}
hipo2(3, 4, angle=pi/2)
hipo(3, 4)
hipo2(3, 4, angle=pi/2+0.3)
hipo2(3, 4, angle=pi/2-0.3)
##-----------------------------------------------------------------------------
## Ângulo com valor default de 90 graus (pi/2 radianos).
hipo3 <- function(a, b, angle=pi/2){
h2 <- a^2+b^2-2*cos(angle)
return(sqrt(h2))
}
hipo3(3, 4)
hipo3(3, 4, angle=pi/2+0.5)
##-----------------------------------------------------------------------------
## Três implementações.
baskara1 <- function(a, b, c){
x1 <- (-b-sqrt(b^2-4*a*c))/(2*a)
x2 <- (-b+sqrt(b^2-4*a*c))/(2*a)
return(c(x1, x2))
}
baskara2 <- function(a, b, c){
s <- sqrt(b^2-4*a*c)
d <- 2*a
x1 <- (-b-s)/d
x2 <- (-b+s)/d
return(c(x1, x2))
}
baskara3 <- function(a, b, c){
x <- (-b+c(-1,1)*sqrt(b^2-4*a*c))/(2*a)
return(x)
}
a <- 2; b <- 1; c <- -3
curve(a*x^2+b*x+c, -3, 3)
abline(h=0, lty=2)
baskara1(a, b, c)
baskara2(a, b, c)
baskara3(a, b, c)
##-----------------------------------------------------------------------------
## Tempo para 50 mil excussões de cada uma.
system.time(replicate(50000, baskara1(a, b, c)))
system.time(replicate(50000, baskara2(a, b, c)))
system.time(replicate(50000, baskara3(a, b, c)))
## As implementação diferentes implicam em custos diferentes.
##-----------------------------------------------------------------------------
## Parábola sem raízes.
a <- 2; b <- 1; c <- 3
curve(a*x^2+b*x+c, -5, 5)
abline(h=0, lty=2)
baskara1(a, b, c)
##-----------------------------------------------------------------------------
## a==0, então não tem curvatura.
a <- 0; b <- 1; c <- 3
curve(a*x^2+b*x+c, -5, 5)
baskara1(a, b, c)
## Como notificar da origem desses erros?
##-----------------------------------------------------------------------------
## Colocando mensagens de erro para testes feitos nos argumentos.
baskara4 <- function(a, b, c){
if (a==0){
stop("O argumento `a` tem que ser diferente de zero.")
}
delta <- b^2-4*a*c
if (delta<=0) {
stop("Função sem raíz pois `delta` é não positivo.")
}
x <- (-b+c(-1,1)*sqrt(delta))/(2*a)
return(x)
}
baskara4(2, 1, -3)
baskara4(0, 1, -3)
baskara4(2, 1, 3)
##-----------------------------------------------------------------------------
## Função que calcula as raízes da equação e ainda os valores no ponto
## crítico (x, y). O ponto crítico sempre existe se a!=0. As raízes
## podem não existir.
baskara5 <- function(a, b, c){
if (a==0){
stop("O argumento `a` tem que ser diferente de zero.")
}
d <- 2*a
delta <- b^2-4*a*c
if (delta<=0){
warning("Função sem raíz pois `delta` é não positivo.")
r <- c(NA, NA)
} else {
r <- (-b+c(-1,1)*sqrt(delta))/d
}
x <- -b/d
y <- a*x^2+b*x+c
L <- list(x=x, y=y, r=r)
class(L) <- "baskara"
return(L)
}
a <- baskara5(2, 1, -3)
str(a)
curve(2*x^2+1*x-3, -2, 2)
abline(v=a$x, h=a$y, col=2, lty=2)
abline(v=a$r, h=0, col=3, lty=2)
plot.baskara <- function(a){
abline(v=a$x, h=a$y, col=2, lty=2)
abline(v=a$r, h=0, col=3, lty=2)
}
curve(2*x^2+1*x-3, -2, 2)
plot(a)
##-----------------------------------------------------------------------------
## Argumentos como vetores, usa posição.
baskara6 <- function(abc){
if (length(abc)==3L){
x <- (-abc[2]+c(-1,1)*sqrt(abc[2]^2-4*abc[1]*abc[3]))/(2*abc[1])
return(x)
} else {
stop("O vetor `abc` deve ter comprimento 3.")
}
}
baskara6(c(2,1,-3))
baskara6(c(2,1,-3,500))
baskara6(c(2,1))
##-----------------------------------------------------------------------------
## Argumentos como vetores nomeados.
baskara7 <- function(abc){
if (length(abc)==3L & all(names(abc)%in%c("a","b","c"))){
x <- (-abc["b"]+c(-1,1)*
sqrt(abc["b"]^2-4*abc["a"]*abc["c"]))/(2*abc["a"])
return(x)
} else {
stop("O vetor `abc` deve ter comprimento 3 e ser nomeado.")
}
}
baskara7(c(a=2, b=1, c=-3))
baskara7(c(b=1, a=2, c=-3))
baskara7(c(b=1, a=2, m=-3))
##-----------------------------------------------------------------------------
## Argumentos como vetores ou lista nomeados.
baskara8 <- function(abc){
if (is.vector(abc)) abc <- as.list(abc)
if (is.list(abc) & length(abc)==3L & all(names(abc)%in%c("a","b","c"))){
x <- with(abc,
(-b+c(-1,1)*sqrt(b^2-4*a*c))/(2*a))
return(x)
} else {
stop("O objeto `abc` deve ser vetor/lista nomeado de 3 elementos.")
}
}
baskara8(c(a=2, b=1, c=-3))
baskara8(list(a=2, b=1, c=-3))
baskara8(data.frame(a=2, b=1, c=-3))print(sessionInfo(), locale=FALSE)## R version 3.1.1 (2014-07-10)
## Platform: i686-pc-linux-gnu (32-bit)
##
## attached base packages:
## [1] methods stats graphics grDevices utils datasets base
##
## other attached packages:
## [1] EnvStats_1.0.3 fBasics_3011.87 timeSeries_3010.97 timeDate_3010.98
## [5] rmarkdown_0.3.3 knitr_1.7
##
## loaded via a namespace (and not attached):
## [1] digest_0.6.4 evaluate_0.5.5 formatR_1.0 htmltools_0.2.6 stringr_0.6.2
## [6] tools_3.1.1 yaml_2.1.13Sys.time()## [1] "2014-11-12 22:04:40 BRST"