Ejercicio 1. (Adaptación ejercicio 15)

Para los datos: 1.0,2.3,4.2,7.1,10.4, utilizar el procedimiento más adecuado para contrastar la hipótesis nula:

  1. Exponencial, \(F_0(x)=1-\exp(-\lambda x), x>0\). Tanto a partir de los valores críticos de las tablas (en papel); como a partir de los valores críticos simulados para distintos \(n\) y \(\alpha\) y concluir a partir de ellos.
# Datos
x<-c(1.0,2.3,4.2,7.1,10.4)

# H0 compuesta
#lillie.test unicamente para normalidad. 

# Para obtener el estadistico podemos. Tambien nos da el pvalor, pero..
emv<-1/mean(x)
ks.test(x,pexp,emv) # Solo puedo usar estadistico. P-valor con la tabla simulada (Exp)
## 
##  Exact one-sample Kolmogorov-Smirnov test
## 
## data:  x
## D = 0.18127, p-value = 0.9864
## alternative hypothesis: two-sided
ks.test(x/mean(x),pexp) # Mismo valor
## 
##  Exact one-sample Kolmogorov-Smirnov test
## 
## data:  x/mean(x)
## D = 0.18127, p-value = 0.9864
## alternative hypothesis: two-sided
#Opcion 1: Obtenido el estadistico podemos ver taba Lilliefors para exponencial para obtener p-valor

#Opcion 2: Obtenido el estadistico podemos comparar con valores criticos simulados para Lilliefors (Exp)
# Tabla de Lilliefors para exponencial

# Modificada para n y alpha
enes<-c(4:12)
alphas<-c(0.95,0.9,0.75,0.5,0.25,0.10,0.075,0.05,0.001)
DnExp<-matrix(0,length(enes),length(alphas))
rownames(DnExp)<-as.character(enes)
colnames(DnExp)<-as.character(alphas)

fila<-0
for(k in enes){
    fila<-fila+1
    columna<-0
    for(j in alphas){
        columna<-columna+1
        matrixExp<-matrix(0,1000,k)
        Dns<-c()
        for(i in 1:1000){
            matrixExp[i,]<-rexp(k,emv) # Si en lugar de emv=0.2, hubiese puesto 5, que ocurriria
            xBar<-mean(matrixExp[i,])
            z<-matrixExp[i,]/xBar
            Dns<-c(Dns,ks.test(z,pexp)$statistic)
            #Dns<-c(Dns,ks.test(x,pexp,1/xBar)$statistic) # O tambien asi
        }
        DnExp[fila,columna]<-quantile(Dns,probs=1-j)
    }
}
DnExp
##         0.95       0.9      0.75       0.5      0.25       0.1     0.075
## 4  0.2022052 0.2314901 0.2681363 0.3201594 0.3824343 0.4508463 0.4740505
## 5  0.1868796 0.2007170 0.2414417 0.2881249 0.3456596 0.4071357 0.4081888
## 6  0.1724985 0.1908508 0.2250361 0.2639116 0.3231668 0.3775028 0.3855225
## 7  0.1601193 0.1748627 0.2011234 0.2447000 0.3042004 0.3518813 0.3594509
## 8  0.1490446 0.1642897 0.1904539 0.2347311 0.2766680 0.3319777 0.3393940
## 9  0.1443641 0.1544561 0.1815257 0.2172664 0.2674339 0.3142067 0.3163489
## 10 0.1359808 0.1481803 0.1710621 0.2119805 0.2475361 0.2987709 0.3085247
## 11 0.1300156 0.1437637 0.1668884 0.1973389 0.2340024 0.2760235 0.2928681
## 12 0.1248593 0.1343412 0.1577496 0.1909576 0.2301706 0.2717661 0.2849452
##         0.05     0.001
## 4  0.4847267 0.6204248
## 5  0.4351619 0.5480261
## 6  0.4092320 0.5481051
## 7  0.3715740 0.5126685
## 8  0.3646249 0.5089865
## 9  0.3514596 0.4343554
## 10 0.3195415 0.4702234
## 11 0.3150971 0.4288154
## 12 0.3024922 0.3888107
  1. Calcular el p-valor por simulación para el test dado en a).
xBar<-mean(x)
tobs<-ks.test(x,pexp,1/xBar)$statistic
DnsH0<-c()
simExp<-matrix(0,1000,length(x))
# Lo repito como si no hubieramos hecho a), pero bastaria con considerar Dns
for(i in 1:1000){
  simExp[i,]<-rexp(length(x),1/xBar)
    DnsH0<-c(DnsH0,ks.test(simExp[i,],pexp,1/simExp[i,])$statistic)
}
pvalExp<-sum(DnsH0>tobs)/1000
pvalExp
## [1] 1
  1. Normal. Tanto desde funciones específicas de R, como a partir de los valores críticos de las tablas (en papel); y también desde de los valores críticos simulados para distintos \(n\) y \(\alpha\) y concluir a partir de ellos.
# H0 compuesta 
m<-mean(x)
s<-sd(x)

# Datos tipificados
z<-(x-m)/s

# Podemos hacerlo varias formas
# 1.- Usando el estadistico y p-valor de lillie.test de la libreria normtest
library(nortest)
lillie.test(x) # El p-valor sirve. Mismo valor
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  x
## D = 0.18356, p-value = 0.8379
lillie.test(z) # El p-valor sirve
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  z
## D = 0.18356, p-value = 0.8379
# 2.- Usando el estadistico de ks.test (para el estadistico) y las tablas Lilliefors (para pvalor)
ks.test(x,pnorm,m,s) # El p-valor no sirve. Mismo valor
## 
##  Exact one-sample Kolmogorov-Smirnov test
## 
## data:  x
## D = 0.18356, p-value = 0.9844
## alternative hypothesis: two-sided
ks.test(z,pnorm) # El p-valor no sirve 
## 
##  Exact one-sample Kolmogorov-Smirnov test
## 
## data:  z
## D = 0.18356, p-value = 0.9844
## alternative hypothesis: two-sided
# 3.- Usando el estadistico de ks.test y las tablas Lilliefors simuladas (Norm)
enes<-c(4:12)
alphas<-c(0.95,0.9,0.85,0.8,0.75,0.5,0.25,0.1,0.05,0.01)
DnNorm<-matrix(0,length(enes),length(alphas))
rownames(DnNorm)<-as.character(enes)
colnames(DnNorm)<-as.character(alphas)
set.seed(1234)
fila<-0
for(k in enes){
    fila<-fila+1
    columna<-0
    for(j in alphas){
        columna<-columna+1
        matrixNorm<-matrix(0,1000,k)
        Dns<-c()
        for(i in 1:1000){
            matrixNorm[i,]<-rnorm(k,m,s)
            xBar<-mean(matrixNorm[i,])
            st<-sd(matrixNorm[i,])
            z<-(matrixNorm[i,]-xBar)/st
            Dns<-c(Dns,ks.test(z,pnorm)$statistic) # Se podria haber usado lillie.test
        }
        DnNorm[fila,columna]<-quantile(Dns,probs=1-j)
    }
}
DnNorm
##         0.95       0.9      0.85       0.8      0.75       0.5      0.25
## 4  0.1669930 0.1866025 0.1955884 0.2115697 0.2192155 0.2579405 0.2943471
## 5  0.1641690 0.1791922 0.1835936 0.1936115 0.2025611 0.2311249 0.2756100
## 6  0.1496122 0.1662484 0.1708479 0.1769048 0.1854476 0.2200505 0.2579886
## 7  0.1410891 0.1522087 0.1624851 0.1666594 0.1746485 0.2083721 0.2413215
## 8  0.1354395 0.1451497 0.1516915 0.1599003 0.1610306 0.1946597 0.2302603
## 9  0.1274311 0.1395766 0.1425364 0.1521915 0.1562692 0.1849371 0.2176267
## 10 0.1182865 0.1318627 0.1398881 0.1440289 0.1485597 0.1780265 0.2054360
## 11 0.1168919 0.1253377 0.1335612 0.1395344 0.1446921 0.1665040 0.1989754
## 12 0.1104904 0.1215888 0.1257444 0.1321678 0.1373488 0.1644175 0.1926626
##          0.1      0.05      0.01
## 4  0.3389818 0.3719604 0.4138511
## 5  0.3176857 0.3411115 0.4011310
## 6  0.3019926 0.3198396 0.3727551
## 7  0.2797531 0.2956719 0.3439593
## 8  0.2671091 0.2799000 0.3325922
## 9  0.2485654 0.2717517 0.3127179
## 10 0.2424919 0.2546090 0.2897199
## 11 0.2283164 0.2465327 0.2875371
## 12 0.2230690 0.2505217 0.2749604
###############################
#        NOTA                 #
###############################
# La libreia nortest permite calcular otros test
# cvm.test() # Cramer von Mises , n>7
# ad.test()  # Anderson Darling , n>7
# sf.test()   # Shapiro Francia, se rechaza para valores pequenos
# shapiro.test() # No esta en la libreria nortest
  1. Calcular por simulación la potencia para la alternativa Exp(7) para el test dado en c) de nivel \(\alpha=0.05\).
vc<-0.344 # A partir de la tabla simulada en c). Tambien 0.344 de la tabla Lilliefors (Norm)
DnsH1<-c()
simExp1<-matrix(0,1000,length(x))
for(i in 1:1000){
  simExp1[i,]<-rexp(length(x),7)
    DnsH1<-c(DnsH1,lillie.test(simExp1[i,])$statistic)
}
potencia<-sum(DnsH1>vc)/1000
potencia
## [1] 0.122
# ¿Como puede conseguirse aumentar potencia? 2 maneras. Disminuir alpha (no recomendado) o aumentar n


###############################
#        NOTA                 #
###############################
# Los valores criticos, pvalores y potencias pueden simularse para cualquier distribucion que conozcamos el estadistico

Ejercicio 2. (Adaptación ejercicio 16)

Diez estudiantes se someten a un test. Los resultados (sobre 100) del test son 95,80,40,52,60,80,82,58,65,50.

  1. Contrastar la hipótesis nula \(H_0:F_0(x)= \left\{ \begin{array}{lcc} 0 & si & x < 0 \\ \\ x^2(3-2x) & si & 0 \leq x \leq 1 \\ \\ 1 & si & x > 1 \end{array} \right.\), considere \(\alpha=0.05.\)
# Necesitamos los datos sobre 1 punto
x<-c(95,80,40,52,60,80,82,58,65,50)/100
n<-length(x)

F0<-function(z){
  return(z^2*(3-2*z))
}

#ordenamos lo datos y calculo F0 para datos ordendos
xo<-sort(x)
F0(xo)
##  [1] 0.352000 0.500000 0.529984 0.618976 0.648000 0.718250 0.896000 0.896000
##  [9] 0.914464 0.992750
plot(ecdf(x)) # Fn
points(xo,F0(xo)) # F0

# ¿Se rechazara H0?

# Estadistico K-S. Calculado directamente con la definicion 
max((1:10)/n-F0(xo))
## [1] 0.00725
max(F0(xo)-(0:9)/n)
## [1] 0.4
Dn<-max(max((1:10)/n-F0(xo)),max(F0(xo)-(0:9)/n))
Dn
## [1] 0.4
# Calculamos el p-valor usando las tablas 
# Entre 0.05 y 0.1

# Estadistico K-S y pvalor calculado con la funcion ks.test
ks.test(x,F0)
## Warning in ks.test.default(x, F0): ties should not be present for the
## Kolmogorov-Smirnov test
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  x
## D = 0.4, p-value = 0.08152
## alternative hypothesis: two-sided
# Tambien podemos usar solo el estadistico y buscar el p-valor en las tablas K-S menos preciso

# Podemos calcular ese pvalor por simulacion. Mandar como ejercicio de entrega

b.Obtener la potencia del test de K-S del apartado a) por simulación frente a la alternativa Beta(2,4).

# Valor critico del test de K-S tablas o `ks.test() o por simuladcion si lo tuvieramos
vc<-0.409 # De las tablas de la distribucion exacta de K-S

DnH1<-c()
for(i in 1:1000){
    sampleBeta<-rbeta(10,2,4)
    DnH1<-c(DnH1,ks.test(sampleBeta,F0)$statistic)
}
sum(DnH1>vc)/1000
## [1] 0.602
  1. Calcular mediante simulación el tamagno muestral necesario para que la potencia sea al menos de 0.95.
# 1.- Lo hacemos por simulacion, tomamos los valores criticos de las tablas
nc<-c(10,20,30,40,50,60,70,80,90,100)
vcs<-c(0.409,0.294,0.242,0.210,1.36/sqrt(nc[5:10]))

potencia<-c()
for(k in 1:length(nc)){
  DnH1<-c()
  for(i in 1:1000){
    sampleBeta<-rbeta(nc[k],2,4)
    DnH1<-c(DnH1,ks.test(sampleBeta,F0)$statistic)
  }
  potencia[k]<-sum(DnH1>vcs[k])/1000
}
potencia
##  [1] 0.611 0.894 0.990 0.997 1.000 1.000 1.000 1.000 1.000 1.000
# Entre n=20 y n=30, puedo afinar un poco mas entre estos dos valores usando los valores criticos de la tabla
nc<-21:29
vcs<-c(0.287,0.281,0.275,0.269,0.264,0.259,0.254,0.250,0.246)

potencia<-c()
for(k in 1:length(nc)){
  DnH1<-c()
  for(i in 1:1000){
    sampleBeta<-rbeta(nc[k],2,4)
    DnH1<-c(DnH1,ks.test(sampleBeta,F0)$statistic)
  }
  potencia[k]<-sum(DnH1>vcs[k])/1000
}
potencia # n=25
## [1] 0.922 0.923 0.938 0.946 0.961 0.964 0.967 0.974 0.975
  1. Obtener la banda de confianza simultanea para la función de distribución desconocida con confianza 0.9 ¿Contiene a la función de distribución bajo la \(H_0\) dada en el apartado a)?
Fn<-ecdf(x)
plot(Fn)
ceros<-rep(0,100)
unos<-rep(1,100)
z<-seq(0,1,length=100)
L<-pmax(ceros,Fn(z)-0.369) # De donde sale este valor
U<-pmin(unos,Fn(z)+0.369)
lines(z,L,type="s",lty=2)
lines(z,U,type="s",lty=2)
lines(xo,F0(xo),type="s",col=2)

Ejercicio 3. (Ejercicio 13).

Se lanzan 5 monedas en 100 ocasiones y se anota el número de caras en cada lanzamiento. Los resultados fueron:

Nº caras 0 1 2 3 4 5
Frecuencias 1 10 20 36 23 10
  1. Estudiar si los datos soportan o no la hipótesis de que todas las monedas son equilibradas.
# H0: b(5,0.5) simple

n<-100
frecO<-c(1,10,20,36,23,10)
frecE<-c(dbinom(0:5,5,0.5)*n) # 5 o mas caras
frecE # <=5
## [1]  3.125 15.625 31.250 31.250 15.625  3.125
# Agrupamos la primera y la ultima clase
n<-100
frecO<-c(11,20,36,33)
frecE<-c(sum(dbinom(0:1,5,0.5))*n,dbinom(2:3,5,0.5)*n,sum(dbinom(4:5,5,0.5))*n)
frecE
## [1] 18.75 31.25 31.25 18.75
# Como en la teoria 
tObs<-sum((frecO-frecE)^2/frecE)
tObs
## [1] 18.80533
pval<-1-pchisq(tObs,3)
pval
## [1] 0.0002999421
# Directamente
chisq.test(x=frecO,p=frecE/n)
## 
##  Chi-squared test for given probabilities
## 
## data:  frecO
## X-squared = 18.805, df = 3, p-value = 0.0002999
  1. Calcular la potencia en la alternativa Bin(5,0.55). ¿Cuánto vale \(n\) para que la potencia sea de 0.7?
p0<-frecE/n
p1<-c(sum(dbinom(0:1,5,0.55)),dbinom(2:3,5,0.55),sum(dbinom(4:5,5,0.55)))
D<-sum((p0-p1)^2/p0)
delta<-n*D
1-pchisq(qchisq(0.95,3),3,delta)
## [1] 0.4270591
# Calculamos n, usando las tablas
deltaNuevo<-8.792
nNuevo<-deltaNuevo/D
ceiling(nNuevo)
## [1] 182
  1. Estudiar si los datos soportan o no la hipótesis de que todas las monedas tienen la misma probabilidad de cara.
# H0: b(5,p) compuesta

# EMV para p en base a los datos agrupados
n<-100
frecO<-c(11,20,36,33)

mlv<-function(p){
  return(-sum(frecO*log(c(sum(dbinom(0:1,5,p)),dbinom(2:3,5,p),sum(dbinom(4:5,5,p))))))
}
minimo<-optim(0.5,mlv,method="Brent",lower=0.01,upper=0.99)
emv<-minimo$par


p0<-c(sum(dbinom(0:1,5,emv)),dbinom(2:3,5,emv),sum(dbinom(4:5,5,emv)))  

# Calculamos el estadistico y el pvalor
tobs<-chisq.test(x=frecO,p=p0)$statistic # Sirve este pvalor...
1-pchisq(tobs,6-1-1-1-1)
## X-squared 
##  0.611787

Ejercicio 4. (Ejercicio 12).

En la generación de 100 observaciones de una distribución Binomial Negativa [BN(2,2/3)], se obtuvieron los resultados siguientes:

Valores 0 1 2 3 4 o mas
Frecuencias 40 24 16 14 6
  1. Contrastar el ajuste a una distribución de Poisson de media 1 (ignoramos binomial negativa).
n<-100
frecO<-c(40,24,16,14,6)
frecE<-c(dpois(0:3,1)*n,(1-ppois(3,1))*n) 
frecE # <=5
## [1] 36.787944 36.787944 18.393972  6.131324  1.898816
chisq.test(x=frecO,p=frecE/n)
## Warning in chisq.test(x = frecO, p = frecE/n): Chi-squared approximation may be
## incorrect
## 
##  Chi-squared test for given probabilities
## 
## data:  frecO
## X-squared = 23.994, df = 4, p-value = 8.011e-05
# Agrupamos la ultima clase n<=5
  1. Repetir el ajuste, agrupando clases hasta que el nº esperado de observaciones en cada clase sea mayor que 5.
n<-100
frecO<-c(40,24,16,20)
frecE<-c(dpois(0:2,1)*n,(1-ppois(2,1))*n) 
frecE 
## [1] 36.78794 36.78794 18.39397  8.03014
chisq.test(x=frecO,p=frecE/n)
## 
##  Chi-squared test for given probabilities
## 
## data:  frecO
## X-squared = 22.88, df = 3, p-value = 4.278e-05
# Agrupamos la ultima clase n<=10
  1. Repetir el ajuste, agrupando clases hasta que el nº esperado de observaciones en cada clase sea mayor que 10.
n<-100
frecO<-c(40,24,36)
frecE<-c(dpois(0:1,1)*n,(1-ppois(1,1))*n)
frecE 
## [1] 36.78794 36.78794 26.42411
chisq.test(x=frecO,p=frecE/n)
## 
##  Chi-squared test for given probabilities
## 
## data:  frecO
## X-squared = 8.1959, df = 2, p-value = 0.01661
  1. Estudiar el efecto que tiene sobre la potencia, en la alternativa BN(2,2/3), considerar unas clases u otras. Tomar el nivel 0.05 para el test.
p0_5clases<-c(dpois(0:3,1),(1-ppois(3,1))) 
p0_4clases<-c(dpois(0:2,1),(1-ppois(2,1))) 
p0_3clases<-c(dpois(0:1,1),(1-ppois(1,1))) 

p1_5clases<-c(dnbinom(0:3,2,2/3),(1-pnbinom(3,2,2/3))) 
p1_4clases<-c(dnbinom(0:2,2,2/3),(1-pnbinom(2,2,2/3))) 
p1_3clases<-c(dnbinom(0:1,2,2/3),(1-pnbinom(1,2,2/3))) 

D_5clases<-sum(((p0_5clases-p1_5clases)^2)/p0_5clases)
D_4clases<-sum(((p0_4clases-p1_4clases)^2)/p0_4clases)
D_3clases<-sum(((p0_3clases-p1_3clases)^2)/p0_3clases)

delta_5clases<-n*D_5clases
delta_4clases<-n*D_4clases
delta_3clases<-n*D_3clases

vc_5clases<-qchisq(0.95,4)
vc_4clases<-qchisq(0.95,3)
vc_3clases<-qchisq(0.95,2)

1-pchisq(vc_5clases,4,delta_5clases)
## [1] 0.5635264
1-pchisq(vc_4clases,3,delta_4clases)
## [1] 0.4296486
1-pchisq(vc_3clases,2,delta_3clases)
## [1] 0.3211214
  1. ¿Cuál debe ser el tamaño muestral para que la potencia sea al menos 0.9, utilizando 5 clases?
deltaNuevo_5clases<-15.405 # Tabla chi2 descentrada
nNuevo_5clases<-deltaNuevo_5clases/D_5clases
ceiling(nNuevo_5clases)
## [1] 210

Ejercicio Propuesto 1. (Adaptación ejercicio 17).

Los diámetros, en cm, de nueve arandelas metálicas, fabricadas por una máquina fueron: 0.962, 1.066, 0.900, 0.846,0.807, 0.797, 0.814, 0.710, 0.676.

  1. ¿Estos datos corresponden a una distribución normal? Calcular el pvalor por simulación.

  2. Calcular por simulación la potencia para la alternativa N(0.8,0.1).

  3. Obtener el tamagno muestral para que la potencia calculada en b) sea de 0.95

Ejercicio Propuesto 2.

Ejercicio 6 de la hoja de ejercicios de TBA.

Ejercicio Propuesto 3.

Obtener por simlación la distribución bajo \(H_0\) del estadístico de Cramer von Mises.