In [1]:
# rm(list=ls())
options(OutDec = ",")
#===============================================================================
silence <- suppressPackageStartupMessages # Para omitir mensagens de alertas
silence(library(psych)) # Para a funcao 'describe'
In [2]:
#===============================================================================
# Seja V o máximo e U o mínimo de uma amostra aleatória simples de N(0,1)
#===============================================================================
k <- 2 #tamanho da amostra
n <- 100000
x <- matrix(rnorm(k*n),n,k)
u <- apply(x,1,"min")
v <- apply(x,1,"max")
r <- v-u
In [3]:
#===============================================================================
# Gráfico de dispersão
#===============================================================================
par(mfrow=c(1,1),lwd=2.0,cex.lab=1.5,cex.axis=1.5,lab=c(7,7,5),
mar=c(5,5,2,2.5),xpd=T,cex.main=2.0,bty="n")
plot(u,v,pch=".",col="red",xlim=c(-6,6),ylim=c(-6,6),axes=F)
axis(side=2,pos=c(0,0))
axis(side=1,pos=c(0,0))
lines(-6:6,-6:6,col="blue",lwd=2,lty=2)
In [4]:
#===============================================================================
# Histograma do mínimo
#===============================================================================
par(mfrow=c(1,1),lwd=2.0,cex.lab=1.5,cex.axis=1.5,lab=c(7,6,5),
mar=c(5,5,2,2.5),xpd=T,cex.main=2.0,bty="n")
hist(u,nclass=50,prob=T,main="",xlab=expression(u),ylab="Densidade",
col="darkorange",xlim=c(-6,6),ylim=c(0,0.6))
In [5]:
#===============================================================================
# Histograma do máximo
#===============================================================================
par(mfrow=c(1,1),lwd=2.0,cex.lab=1.5,cex.axis=1.5,lab=c(7,6,5),
mar=c(5,5,2,2.5),xpd=T,cex.main=2.0,bty="n")
hist(v,nclass=50,prob=T,main="",xlab=expression(v),ylab="Densidade",
col="darkorange",xlim=c(-6,6),ylim=c(0,0.6))
In [6]:
#===============================================================================
# Histograma da amplitude
#===============================================================================
par(mfrow=c(1,1),lwd=2.0,cex.lab=1.5,cex.axis=1.5,lab=c(7,6,5),
mar=c(5,5,2,2.5),xpd=T,cex.main=2.0,bty="n")
hist(r,nclass=50,prob=T,main="",xlab=expression(r),ylab="Densidade",
col="darkblue",xlim=c(0,8),ylim=c(0,0.6))
In [7]:
#===============================================================================
# Estatísticas descritivas da amplitude
#===============================================================================
options(digits=5)
rdescribe <- describe(r)
rmedia <- rdescribe[,3]
rdesvpad <- rdescribe[,4]
rmediana <- rdescribe[,5]
rmin <- rdescribe[,8]
rmax <- rdescribe[,9]
rrange <- rdescribe[,10]
rassimetria <- rdescribe[,11]
rexccurtose <- rdescribe[,12]
STATS <- cbind(rmedia,rdesvpad,rmediana,rmin,rmax,rrange,
rassimetria,rexccurtose)
print(STATS)
rmedia rdesvpad rmediana rmin rmax rrange rassimetria rexccurtose [1,] 1,1277 0,85339 0,95133 0,0001017 5,8608 5,8607 0,99762 0,86183
In [8]:
#===============================================================================
# Início da avaliação das constantes via aproximação de Riemann
#===============================================================================
n <- seq(2,40,1) # sequência de tamanhos da amostras
m <- length(n) # número de tamnhos de amostras
In [9]:
#===============================================================================
# Analisando a distribuição do mínimo
#===============================================================================
EU <- rep(NA,m) # valor esperado de U
EU2 <- rep(NA,m) # valor esperado de U^2
M <- 10000000 # tamanho da grade na aproximação de Riemann
lbu <- -9
ubu <- 9
useq <- seq(lbu,ubu,length=M)
incu <- useq[2]-useq[1]
for(j in 1:m){
aux <- useq*n[j]*dnorm(useq)*((1-pnorm(useq))^(n[j]-1))
EU[j] <- sum(aux)*incu
EU2[j] <- sum(useq*aux)*incu
}
VARU <- EU2 - EU^2
print(cbind(EU,VARU))
rm(list=c("aux","useq"))
EU VARU [1,] -0,56419 0,68169 [2,] -0,84628 0,55947 [3,] -1,02938 0,49172 [4,] -1,16296 0,44753 [5,] -1,26721 0,41593 [6,] -1,35218 0,39192 [7,] -1,42360 0,37290 [8,] -1,48501 0,35735 [9,] -1,53875 0,34434 [10,] -1,58644 0,33325 [11,] -1,62923 0,32364 [12,] -1,66799 0,31521 [13,] -1,70338 0,30773 [14,] -1,73591 0,30104 [15,] -1,76599 0,29501 [16,] -1,79394 0,28953 [17,] -1,82003 0,28453 [18,] -1,84448 0,27994 [19,] -1,86748 0,27570 [20,] -1,88917 0,27177 [21,] -1,90969 0,26811 [22,] -1,92916 0,26470 [23,] -1,94767 0,26151 [24,] -1,96531 0,25851 [25,] -1,98216 0,25569 [26,] -1,99827 0,25302 [27,] -2,01371 0,25050 [28,] -2,02852 0,24811 [29,] -2,04276 0,24584 [30,] -2,05646 0,24368 [31,] -2,06967 0,24162 [32,] -2,08241 0,23965 [33,] -2,09471 0,23777 [34,] -2,10661 0,23597 [35,] -2,11812 0,23425 [36,] -2,12928 0,23259 [37,] -2,14009 0,23100 [38,] -2,15059 0,22947 [39,] -2,16078 0,22799
In [10]:
#===============================================================================
# Analisando a distribuição do máximo
#===============================================================================
EV <- rep(NA,m) # valor esperado de V
EV2 <- rep(NA,m) # valor esperado de V^2
N <- 10000000 # tamanho da grade na aproximação de Riemann
lbv <- -9
ubv <- 9
vseq <- seq(lbv,ubv,length=N)
incv <- vseq[2]-vseq[1]
for(j in 1:m){
aux <- vseq*n[j]*dnorm(vseq)*(pnorm(vseq)^(n[j]-1))
EV[j] <- sum(aux)*incv
EV2[j] <- sum(vseq*aux)*incv
}
VARV <- EV2 - EV^2
print(cbind(EV,VARV))
rm(list=c("aux","vseq"))
EV VARV [1,] 0,56419 0,68169 [2,] 0,84628 0,55947 [3,] 1,02938 0,49172 [4,] 1,16296 0,44753 [5,] 1,26721 0,41593 [6,] 1,35218 0,39192 [7,] 1,42360 0,37290 [8,] 1,48501 0,35735 [9,] 1,53875 0,34434 [10,] 1,58644 0,33325 [11,] 1,62923 0,32364 [12,] 1,66799 0,31521 [13,] 1,70338 0,30773 [14,] 1,73591 0,30104 [15,] 1,76599 0,29501 [16,] 1,79394 0,28953 [17,] 1,82003 0,28453 [18,] 1,84448 0,27994 [19,] 1,86748 0,27570 [20,] 1,88917 0,27177 [21,] 1,90969 0,26811 [22,] 1,92916 0,26470 [23,] 1,94767 0,26151 [24,] 1,96531 0,25851 [25,] 1,98216 0,25569 [26,] 1,99827 0,25302 [27,] 2,01371 0,25050 [28,] 2,02852 0,24811 [29,] 2,04276 0,24584 [30,] 2,05646 0,24368 [31,] 2,06967 0,24162 [32,] 2,08241 0,23965 [33,] 2,09471 0,23777 [34,] 2,10661 0,23597 [35,] 2,11812 0,23425 [36,] 2,12928 0,23259 [37,] 2,14009 0,23100 [38,] 2,15059 0,22947 [39,] 2,16078 0,22799
In [11]:
#===============================================================================
# Como Z ~ N(0,1), compara-se os resultados obtidos para mínimo e máximo
# Note que E(V) = -E(U) e Var(V) = Var(U).
#===============================================================================
compara <- cbind(EU, VARU, EV, VARV)
print(compara)
EU VARU EV VARV [1,] -0,56419 0,68169 0,56419 0,68169 [2,] -0,84628 0,55947 0,84628 0,55947 [3,] -1,02938 0,49172 1,02938 0,49172 [4,] -1,16296 0,44753 1,16296 0,44753 [5,] -1,26721 0,41593 1,26721 0,41593 [6,] -1,35218 0,39192 1,35218 0,39192 [7,] -1,42360 0,37290 1,42360 0,37290 [8,] -1,48501 0,35735 1,48501 0,35735 [9,] -1,53875 0,34434 1,53875 0,34434 [10,] -1,58644 0,33325 1,58644 0,33325 [11,] -1,62923 0,32364 1,62923 0,32364 [12,] -1,66799 0,31521 1,66799 0,31521 [13,] -1,70338 0,30773 1,70338 0,30773 [14,] -1,73591 0,30104 1,73591 0,30104 [15,] -1,76599 0,29501 1,76599 0,29501 [16,] -1,79394 0,28953 1,79394 0,28953 [17,] -1,82003 0,28453 1,82003 0,28453 [18,] -1,84448 0,27994 1,84448 0,27994 [19,] -1,86748 0,27570 1,86748 0,27570 [20,] -1,88917 0,27177 1,88917 0,27177 [21,] -1,90969 0,26811 1,90969 0,26811 [22,] -1,92916 0,26470 1,92916 0,26470 [23,] -1,94767 0,26151 1,94767 0,26151 [24,] -1,96531 0,25851 1,96531 0,25851 [25,] -1,98216 0,25569 1,98216 0,25569 [26,] -1,99827 0,25302 1,99827 0,25302 [27,] -2,01371 0,25050 2,01371 0,25050 [28,] -2,02852 0,24811 2,02852 0,24811 [29,] -2,04276 0,24584 2,04276 0,24584 [30,] -2,05646 0,24368 2,05646 0,24368 [31,] -2,06967 0,24162 2,06967 0,24162 [32,] -2,08241 0,23965 2,08241 0,23965 [33,] -2,09471 0,23777 2,09471 0,23777 [34,] -2,10661 0,23597 2,10661 0,23597 [35,] -2,11812 0,23425 2,11812 0,23425 [36,] -2,12928 0,23259 2,12928 0,23259 [37,] -2,14009 0,23100 2,14009 0,23100 [38,] -2,15059 0,22947 2,15059 0,22947 [39,] -2,16078 0,22799 2,16078 0,22799
In [12]:
#===============================================================================
# Cálculo da constante d_2 = E(W) = E(V) - E(U) = E(max X_i) - E(min X_i)
#===============================================================================
d2 <- EV - EU # Valor esperado de W = Z_(n) - Z_(1)
print(d2)
[1] 1,1284 1,6926 2,0588 2,3259 2,5344 2,7044 2,8472 2,9700 3,0775 3,1729 [11] 3,2585 3,3360 3,4068 3,4718 3,5320 3,5879 3,6401 3,6890 3,7350 3,7783 [21] 3,8194 3,8583 3,8953 3,9306 3,9643 3,9965 4,0274 4,0570 4,0855 4,1129 [31] 4,1393 4,1648 4,1894 4,2132 4,2362 4,2586 4,2802 4,3012 4,3216
In [13]:
#===============================================================================
# Note que podemos calcular d_2 = 2*E(V) ou d_2 = -2*E(U)
#===============================================================================
print(cbind(d2,2*EV,-2*EU))
d2 [1,] 1,1284 1,1284 1,1284 [2,] 1,6926 1,6926 1,6926 [3,] 2,0588 2,0588 2,0588 [4,] 2,3259 2,3259 2,3259 [5,] 2,5344 2,5344 2,5344 [6,] 2,7044 2,7044 2,7044 [7,] 2,8472 2,8472 2,8472 [8,] 2,9700 2,9700 2,9700 [9,] 3,0775 3,0775 3,0775 [10,] 3,1729 3,1729 3,1729 [11,] 3,2585 3,2585 3,2585 [12,] 3,3360 3,3360 3,3360 [13,] 3,4068 3,4068 3,4068 [14,] 3,4718 3,4718 3,4718 [15,] 3,5320 3,5320 3,5320 [16,] 3,5879 3,5879 3,5879 [17,] 3,6401 3,6401 3,6401 [18,] 3,6890 3,6890 3,6890 [19,] 3,7350 3,7350 3,7350 [20,] 3,7783 3,7783 3,7783 [21,] 3,8194 3,8194 3,8194 [22,] 3,8583 3,8583 3,8583 [23,] 3,8953 3,8953 3,8953 [24,] 3,9306 3,9306 3,9306 [25,] 3,9643 3,9643 3,9643 [26,] 3,9965 3,9965 3,9965 [27,] 4,0274 4,0274 4,0274 [28,] 4,0570 4,0570 4,0570 [29,] 4,0855 4,0855 4,0855 [30,] 4,1129 4,1129 4,1129 [31,] 4,1393 4,1393 4,1393 [32,] 4,1648 4,1648 4,1648 [33,] 4,1894 4,1894 4,1894 [34,] 4,2132 4,2132 4,2132 [35,] 4,2362 4,2362 4,2362 [36,] 4,2586 4,2586 4,2586 [37,] 4,2802 4,2802 4,2802 [38,] 4,3012 4,3012 4,3012 [39,] 4,3216 4,3216 4,3216
In [14]:
#===============================================================================
# Cálculo da variância de W: Var(W). Na verdade, é o cálculo de E(W^2).
#===============================================================================
inc <- 0.01
lb <- -7
ub <- 7
useq <- seq(lb,ub,inc) # mínimo
vseq <- useq # máximo
nuv <- length(vseq)
EW2 <- rep(NA,m)
#-------------------------------------------------------------------------------
# Note que a integração deveria ser em u e w, mas estamos fazendo em u e v
# porque consideramos a transformada v = u + w. Lembre-se que u <= v.
#-------------------------------------------------------------------------------
for(j in 1:m){
aux <- 0
# Avaliação da função f(u,v) em uma grade de valores
for(i in 1:nuv){
for(k in 1:i){
aux <- aux + n[j]*(n[j]-1)*((vseq[i]-useq[k])^2)*
dnorm(useq[i])*dnorm(vseq[k])*
((pnorm(vseq[i])-pnorm(useq[k]))^(n[j]-2))
}
}
EW2[j] <- sum(aux)*inc*inc
}
In [15]:
#===============================================================================
# Cálculo da constante d_3 = DP(W)
#===============================================================================
d3 <- sqrt(EW2 - d2^2)
print(d3)
[1] 0,85250 0,88837 0,87981 0,86408 0,84804 0,83321 0,81983 0,80783 0,79705 [10] 0,78731 0,77848 0,77042 0,76302 0,75621 0,74991 0,74405 0,73859 0,73348 [19] 0,72869 0,72417 0,71991 0,71589 0,71207 0,70844 0,70499 0,70170 0,69855 [28] 0,69555 0,69267 0,68990 0,68725 0,68470 0,68224 0,67987 0,67759 0,67538 [37] 0,67325 0,67118 0,66919
In [16]:
#===============================================================================
# Cálculo da constante c_4
#
# Evite o cálculo: c4 <- sqrt(2/(n-1))*gamma(n/2) /gamma((n-1)/2), pois pode
# ocorrer problemas numéricos.
#
# A versao abaixo é numericamente mais estável e precisa.
#===============================================================================
c4 <- exp(0.5*(log(2)-log(n-1))+lgamma(0.5*n)-lgamma(0.5*(n-1)))
In [17]:
#===============================================================================
# Combinando os resultados das constantes
#===============================================================================
constantes <- cbind(n,d2,d3,c4)
colnames(constantes) <- c("n","d2(n)","d3(n)","c4(n)")
print(constantes)
n d2(n) d3(n) c4(n) [1,] 2 1,1284 0,85250 0,79788 [2,] 3 1,6926 0,88837 0,88623 [3,] 4 2,0588 0,87981 0,92132 [4,] 5 2,3259 0,86408 0,93999 [5,] 6 2,5344 0,84804 0,95153 [6,] 7 2,7044 0,83321 0,95937 [7,] 8 2,8472 0,81983 0,96503 [8,] 9 2,9700 0,80783 0,96931 [9,] 10 3,0775 0,79705 0,97266 [10,] 11 3,1729 0,78731 0,97535 [11,] 12 3,2585 0,77848 0,97756 [12,] 13 3,3360 0,77042 0,97941 [13,] 14 3,4068 0,76302 0,98097 [14,] 15 3,4718 0,75621 0,98232 [15,] 16 3,5320 0,74991 0,98348 [16,] 17 3,5879 0,74405 0,98451 [17,] 18 3,6401 0,73859 0,98541 [18,] 19 3,6890 0,73348 0,98621 [19,] 20 3,7350 0,72869 0,98693 [20,] 21 3,7783 0,72417 0,98758 [21,] 22 3,8194 0,71991 0,98817 [22,] 23 3,8583 0,71589 0,98870 [23,] 24 3,8953 0,71207 0,98919 [24,] 25 3,9306 0,70844 0,98964 [25,] 26 3,9643 0,70499 0,99005 [26,] 27 3,9965 0,70170 0,99043 [27,] 28 4,0274 0,69855 0,99079 [28,] 29 4,0570 0,69555 0,99111 [29,] 30 4,0855 0,69267 0,99142 [30,] 31 4,1129 0,68990 0,99170 [31,] 32 4,1393 0,68725 0,99197 [32,] 33 4,1648 0,68470 0,99222 [33,] 34 4,1894 0,68224 0,99245 [34,] 35 4,2132 0,67987 0,99268 [35,] 36 4,2362 0,67759 0,99288 [36,] 37 4,2586 0,67538 0,99308 [37,] 38 4,2802 0,67325 0,99327 [38,] 39 4,3012 0,67118 0,99344 [39,] 40 4,3216 0,66919 0,99361
In [18]:
#===============================================================================
# Fim
#===============================================================================