Universidade Federal do Paraná
Curso de Estatística
CE 089 -
Estatística Computacional II - 2014/2
Prof. Dr. Walmes Marques Zeviani
Tabela de conteúdo
##-----------------------------------------------------------------------------
## Lançamento de uma moeda.
## x <- scan()
## dput(x)
## Laçamento da moeda pelo Leonardo.
x <- c(1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 1, 0, 0, 0, 0,
0, 1, 1, 1, 0, 1, 1, 0, 1, 1, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 1,
1, 0, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 1, 1, 1, 0,
0, 1, 0, 1, 1, 0, 0, 0)
## Proporção amostral.
sum(x)/length(x)
## [1] 0.4571
## Número de trocas de face.
o <- sum(abs(diff(x))); o
## [1] 27
## Laçamento de uma moeda verdadeira (p=0.5, lanç indep = sem memória)
sample(0:1, 70, replace=TRUE)
## [1] 1 1 0 1 1 0 0 0 1 0 0 0 0 1 1 1 0 0 0 1 1 0 0 0 0 1 0 1 0 1 1 1 1 0 0 0 1 1 1 0 1 1 1
## [44] 0 0 1 1 0 0 1 0 1 0 1 0 1 1 1 0 1 0 1 1 1 0 0 0 0 0 0
sum(abs(diff(sample(0:1, 70, replace=TRUE))))
## [1] 30
## Função que lança uma moeda verdadeira e retorna o número de
## trocas. Ela reproduz o experimento sob a hipótese nula, ou seja, com
## p=0.5 e lançamentos independentes (sem memória).
moeda <- function(n){
## sum(abs(diff(rbinom(n, 1, 0.5))))
sum(abs(diff(sample(0:1, n, replace=TRUE))))
}
moeda(70)
## [1] 39
## Faz várias execuções do experimento aleatório.
r <- replicate(50000, moeda(70))
## A distribuição amostral da estatística número de trocas.
hist(r, breaks=seq(min(r), max(r)+1, by=1)-0.5, prob=TRUE,
xlab="Número de trocas em 70 lançamentos",
ylab="Densidade", main=NULL)
abline(v=o, col=2)
text(x=o, y=par()$usr[4], label="Estatística observada", srt=90,
adj=c(1.25,-0.25))
plot(ecdf(r))
abline(v=o, col=2)
text(x=o, y=par()$usr[4], label="Estatística observada", srt=90,
adj=c(1.25,-0.25))
## Como a v.a. é discreta.
sum(r<=o)/length(r)
## [1] 0.04532
sum(r<o)/length(r)
## [1] 0.02628
## A distribuição do número de trocas sob H0 é uma Binomial.
i <- 0:69
p <- pbinom(i, size=69, p=0.5)
plot(ecdf(r), verticals=TRUE, cex=NA, main=NULL,
xlab="Número de trocas em 70 lançamentos",
ylab="Probabilidade acumulada")
lines(p~i, type="s", col=2)
legend("right",
legend=c("Distribuição empírica", "Distribuição teórica"),
col=1:2, lty=1, bty="n")
##-----------------------------------------------------------------------------
## Teste para média de duas populações.
require(lattice); require(latticeExtra)
## Loading required package: lattice
## Loading required package: latticeExtra
## Loading required package: RColorBrewer
## as.factor(rep(1:2, each=30))
## factor(rep(1:2, each=30), labels=c("N","S"))
da <- data.frame(local=gl(2, 30, labels=c("N","S")))
da$bico <- with(da, rnorm(length(local),
mean=49+2*as.integer(local),
sd=2))
head(da)
## local bico
## 1 N 48.28
## 2 N 51.09
## 3 N 49.43
## 4 N 49.10
## 5 N 50.63
## 6 N 50.19
densityplot(~bico, groups=local, data=da,
auto.key=list(title="Local da ilha", cex.title=1.1),
ylab="Densidade", xlab="Comprimento do bico das aves")
##-----------------------------------------------------------------------------
## Média amostral observada.
a <- with(da, tapply(bico, local, mean)); a
## N S
## 50.72 52.88
d <- diff(a); d
## S
## 2.158
## Aleatorizando os valores já que sob a hipótese nula não há diferença
## de médias entre as população o que implica dizer que os valores são
## todos provenientes da mesma população.
with(da, tapply(sample(bico, replace=TRUE), local, mean))
## N S
## 52.27 51.63
## Distribuição da diferença de médias.
r <- replicate(50000, diff(with(da, tapply(sample(bico), local, mean))))
head(r)
## S S S S S S
## -0.02463 -0.12364 -0.07986 0.11643 -0.43638 -0.04584
par(mfrow=c(1,2))
plot(density(r))
abline(v=d)
plot(ecdf(r))
abline(v=d); layout(1)
## Cálculo do p-valor empírico.
sum(abs(r)>=d)/length(r)
## [1] 8e-05
##-----------------------------------------------------------------------------
## Teste Monte Carlo para a hipótese de distribuição aleatória de
## pontos.
## Pontos aleatórios dentro de um quadrado unitário (as posições são
## independentes dos demais pontos).
x <- runif(50)
y <- runif(50)
plot(x, y, asp=1, axes=FALSE, ann=FALSE)
rect(0, 0, 1, 1, lty=2)
d <- dist(cbind(x,y))
d <- as.vector(d)
c(length(d), 100*99/2)
## [1] 1225 4950
## Menor distância entre dois pontos.
min(d)
## [1] 0.007768
## Distribução acumulada das distência entre pontos.
plot(ecdf(d), xlab="Distância", ylab="Proporção de pontos")
##-----------------------------------------------------------------------------
## Distribuições de eventos espaciais não independentes.
require(spatstat)
## Loading required package: spatstat
##
## spatstat 1.38-1 (nickname: 'Le Hardi')
## For an introduction to spatstat, type 'beginner'
##
## Attaching package: 'spatstat'
##
## The following object is masked from 'package:lattice':
##
## panel.histogram
## Distribuição regular.
xy <- rSSI(n=50, r=0.1)
## plot(xy)
plot(xy$x, xy$y, asp=1, axes=FALSE, ann=FALSE)
rect(0, 0, 1, 1, lty=2)
## Distância para pontos com distribuição regular.
drg <- as.vector(dist(cbind(xy$x, xy$y)))
## Menor distância entre dois pontos.
min(drg)
## [1] 0.1
## Distribução acumulada das distência entre pontos.
plot(ecdf(drg), xlab="Distância", ylab="Proporção de pontos")
## Distribuição agregada.
zk <- rGaussPoisson(40, 0.1, 0.98)
plot(zk$x, zk$y, asp=1, axes=FALSE, ann=FALSE)
rect(0, 0, 1, 1, lty=2)
## Distância para pontos com distribuição regular.
dag <- as.vector(dist(cbind(zk$x, zk$y)))
## Menor distância entre dois pontos.
min(dag)
## [1] 0.009561
## Distribução acumulada das distência entre pontos.
plot(ecdf(dag), xlab="Distância", ylab="Proporção de pontos")
## Todos os gráficos juntos.
plot(ecdf(d), xlab="Distância", ylab="Proporção de pontos")
lines(ecdf(drg), col=2)
lines(ecdf(dag), col=4)
