Taller

Ejercicio 5.6

En el tema 2 se propuso el análisis de la aleatoriedad de un generador de números pseudo-aleatorios mediante el test de rachas, que se podría implementar repetidamente. Sin embargo, la aproximación asintótica empleada por la rutina runs.test de la librería tseries no es adecuada para tamaños muéstrales pequeños \((n<40)\) y sería preferible utilizar la distribución exacta (o por lo menos utilizar una corrección por continuidad).

a. Analizar el comportamiento del contraste empleando repetidamente el test de rachas, considerando 500 pruebas con muestras de tamaño 10 de una Bernoulli(0.5). ¿Se observa algo extraño?

## Warning: package 'tseries' was built under R version 3.5.3

nsim <- 500
n <- 10
p <- 0.5
estadistico <- numeric(nsim)
pvalor <- numeric(nsim)
for(isim in 1:nsim) {
  v <- rbinom(n,1,p) 
  u <- factor(sign(v))
  tmp <- runs.test(u,alternative="two.sided")
  estadistico[isim] <- tmp$statistic
  pvalor[isim] <- tmp$p.value
}

hist(estadistico,lwd=2,border='red',main='Estadistico')

plot(density(estadistico),lwd=2,col='blue',main = 'Curva de Densidad')

# Se observan dos picos, podemos asumir que la distribución es bimodal. 

b Realiza un programa que permita aproximar por simulación la función de masa de probabilidad del estadístico número de rachas (a partir de valores de una Bernoulli(0.5). Representarla gráficamente y compararla con la densidad normal. Obtener los puntos críticos para contrastar la hipótesis nula de aleatoriedad para \(α=0.01\), \(0.05\) y \(0.1\). ¿Es esta dístribución adecuada para el contraste de aleatoriedad de variables continuas?¿Cual debería ser la probabilidad de obtener una única racha al aplicar el test a una variable continua?

media <- mean(estadistico)
sda <- sd(estadistico)
plot(1:nsim, cumsum(estadistico)/(1:nsim),col='orange',lwd=2, type="l", ylab="Media muestral", 
     xlab="Nº de simulaciones")

plot(density(estadistico),col='red',lwd=2)
curve(dnorm(x),col='blue',lwd=2,add=TRUE)

#Converge una Distribucion N(0,1).

## Media = 0.03997718
##  Desviación Estandar 1.028053

#Puntos Criticos.
alfa <- 0.05
z <- qnorm(1 - alfa/2)
Intervalo1 <- c(media - z*1/sqrt(n),media + z*1/sqrt(n))

alfa <- 0.01
z <- qnorm(1 - alfa/2)
Intervalo2 <- c(media - z*1/sqrt(n),media + z*1/sqrt(n))

alfa <- 0.1
z <- qnorm(1 - alfa/2)
Intervalo3 <- c(media - z*1/sqrt(n),media + z*1/sqrt(n))

{
  cat("Alfa 0.05 =", Intervalo1)
  cat("\n Alfa 0.01 =", Intervalo2)
  cat("\n Alfa 0.1 =", Intervalo3)
}

## Alfa 0.05 = -0.5798178 0.6597722
##  Alfa 0.01 = -0.7745716 0.8545259
##  Alfa 0.1 = -0.4801712 0.5601256

# Es adecuada ya que el estadistico converje en distribucion a una N(0,1)

# La probabilidad debería ser 1. 

c Diseñar una rutina que permita realizar el contraste de aleatoriedad de una variable continua aproximando el p-valor por simulación. Asumir que la distribución del estadístico puede ser asimétrica.

estadistico <- numeric(nsim)
pvalor <- numeric(nsim)


for(isim in 1:nsim) {
  u <- rbinom(n,1,p) 
  u <- factor(sign(u))
  tmp <- runs.test(u,alternative="two.sided")
  estadistico[isim] <- tmp$statistic
  pvalor[isim] <- tmp$p.value
}

mean(pvalor)

## [1] 0.5141267

plot(1:nsim, cumsum(pvalor)/(1:nsim),col='orange',lwd=2, type="l", ylab="Media muestral", 
     xlab="Nº de simulaciones")

d Diseñar una rutina que permita realizar el contraste de aleatoriedad de una variable continua aproximando el p-valor mediante bootstrap.

boot.f <- function(data, indices){
  mean(data[indices])
}

nboot <- 1000
set.seed(1)
stat.boot <- boot(pvalor, boot.f, nboot)
stat.boot

## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = pvalor, statistic = boot.f, R = nboot)
## 
## 
## Bootstrap Statistics :
##      original        bias    std. error
## t1* 0.5141267 -0.0002514539  0.01381938

names(stat.boot)

##  [1] "t0"        "t"         "R"         "data"      "seed"     
##  [6] "statistic" "sim"       "call"      "stype"     "strata"   
## [11] "weights"

hist(stat.boot$t, freq=FALSE)

plot(stat.boot)