Práctica 9 (2016)

Pràctica Catalán
Universidad Universidad Autónoma de Barcelona (UAB)
Grado Estadística Aplicada - 3º curso
Asignatura Simulació, Remostreig i Aplicacions
Año del apunte 2016
Páginas 10
Fecha de subida 30/04/2016
Descargas 6
Subido por

Vista previa del texto

Pra´ctica 9: Bootstrap no param´etrico (II) Anna Olmo, 1363013 April 30, 2016 1 Ejercicio 1: Intervalos de confianza Bootstrap (II) Sheehan et al (1991) estudiaron si la desaparici´on de los dinosaurios se produjo gradualmente o fue un hecho espontaneo al final de la era cret´ acica. Para hacerlo, utilizaron el n´ umero de dinosaurios registrados en Dakota del Norte y Montana para diferentes especies de dinosaurios. Los datos recogen esta informaci´ on: Ceratops Hadro Pachycephalo Tyranno Ornithomo Sauronith Inferior Mitj` a Superior 50 53 19 29 51 7 2 19 1 1 3 2 4 8 1 1 6 3 Ceratops Hadro Pachycephalo Tyranno Ornithomo Sauronith Inferior 0.57471264 0.33333333 0.02298851 0.01149425 0.04597701 0.01149425 Mitj` a 0.37857143 0.36428571 0.13571429 0.02142857 0.05714286 0.04285714 Superior 0.57575758 0.21212121 0.03030303 0.06060606 0.03030303 0.09090909 Para medir la diferencia entre los recuentos se puede utilizar la entrop´ıa o diversidad de Havrda-Charvat de orden α que se define por la expresi´ on: 1 [1 − α−1 donde pi son las correspondientes probabilidades.
HCα = pi α ] > HC<-function(alpha,num){ + ( 1/(alpha-1) ) * ( 1 - sum(dinos.prob[,num]^alpha)) + } Su desviaci´ on est´ andar tiene la siguiente expresi´on: seHC = 1 α nα−1 pi 2α−1 − ( pi )2 > HC.se<-function(alpha,num){ + ( 1/sqrt( nrow(dinos.prob) ) ) * abs( alpha/(alpha-1) ) * + sqrt( sum( dinos.prob[,num]^(2*alpha-1) ) + sum(dinos.prob[,num]^alpha)^2 ) + } A continuaci´ on, da respuesta a las siguientes cuestiones: 1 1.1 Calcula las biodiversidades de Havrda-Charvat para los datos del cret´ acico inferior, medio y superior, y sus errores est´ andares utilizando los datos, y α = 0.5, 2, 4.
• α = 0.5 : > > > > > > inferior0.5<-HC(0.5,1) inferior0.5.se<-HC.se(0.5,1) medio0.5<-HC(0.5,2) medio0.5.se<-HC.se(0.5,2) superior0.5<-HC(0.5,3) superior0.5.se<-HC.se(0.5,3) HC se inferior 1.831826 0.6230674 medio 2.359374 0.4562463 superior 2.230405 0.5043011 •α=2: > > > > > > inferior2<-HC(2,1) inferior2.se<-HC.se(2,1) medio2<-HC(2,2) medio2.se<-HC.se(2,2) superior2<-HC(2,3) superior2.se<-HC.se(2,3) HC se inferior 0.5556877 0.1403815 medio 0.7000000 0.1012339 superior 0.6097337 0.1809752 •α=4: > > > > > > inferior4<-HC(4,1) inferior4.se<-HC.se(4,1) medio4<-HC(4,2) medio4.se<-HC.se(4,2) superior4<-HC(4,3) superior4.se<-HC.se(4,3) HC se inferior 0.2928517 0.04360488 medio 0.3204989 0.01197583 superior 0.2960006 0.05003587 1.2 Haz un boostrap para la biodiversidad de Havrda-Charvat utilizando los datos proporcionados. Calcula los intervalos de confianza bootstrap. Eval´ ua la biodiversidad y su correspondiente intervalo de confianza al 95% en los diferentes per´ıodos del cret´ acido para α = 0.5, 2, 4 (no uses boot).
> boots<-function(alpha,num){ + datos<-sample(dinos.prob[,num],replace=TRUE) 2 + HC<-( 1/(alpha-1) ) * ( 1 - sum(datos^alpha)) + HC.se<-( 1/sqrt(length(datos)) ) * abs( alpha/(alpha-1) ) * + sqrt( sum( datos^(2*alpha-1) ) - sum(datos^alpha)^2 ) + c(HC,HC.se) + } • α = 0.5 : > > + + > > + + > > + + rep.inferior0.5<-t(replicate(1000,boots(0.5,1))) int.inferior0.5<-c(mean(rep.inferior0.5[,1]), mean(rep.inferior0.5[,1])-1.96*mean(na.omit(rep.inferior0.5[,2])), mean(rep.inferior0.5[,1])+1.96*mean(na.omit(rep.inferior0.5[,2]))) rep.medio0.5<-t(replicate(1000,boots(0.5,2))) int.medio0.5<-c(mean(rep.medio0.5[,1]), mean(rep.medio0.5[,1])-1.96*mean(na.omit(rep.medio0.5[,2])), mean(rep.medio0.5[,1])+1.96*mean(na.omit(rep.medio0.5[,2]))) rep.superior0.5<-t(replicate(1000,boots(0.5,3))) int.superior0.5<-c(mean(rep.superior0.5[,1]), mean(rep.superior0.5[,1])-1.96*mean(na.omit(rep.superior0.5[,2])), mean(rep.superior0.5[,1])+1.96*mean(na.omit(rep.superior0.5[,2]))) HC inf sup inferior 1.853412 0.5240025 3.182822 medio 2.358490 1.2878205 3.429159 superior 2.261775 1.1136376 3.409912 •α=2: > > + + > > + + > > + + rep.inferior2<-t(replicate(1000,boots(2,1))) int.inferior2<-c(mean(rep.inferior2[,1]), mean(rep.inferior2[,1])-1.96*mean(na.omit(rep.inferior2[,2])), mean(rep.inferior2[,1])+1.96*mean(na.omit(rep.inferior2[,2]))) rep.medio2<-t(replicate(1000,boots(2,2))) int.medio2<-c(mean(rep.medio2[,1]), mean(rep.medio2[,1])-1.96*mean(na.omit(rep.medio2[,2])), mean(rep.medio2[,1])+1.96*mean(na.omit(rep.medio2[,2]))) rep.superior2<-t(replicate(1000,boots(2,3))) int.superior2<-c(mean(rep.superior2[,1]), mean(rep.superior2[,1])-1.96*mean(na.omit(rep.superior2[,2])), mean(rep.superior2[,1])+1.96*mean(na.omit(rep.superior2[,2]))) HC inf sup inferior 0.5650242 0.2935368 0.8365116 medio 0.7017929 0.4938880 0.9096978 superior 0.6131001 0.3724844 0.8537158 •α=4: > rep.inferior4<-t(replicate(1000,boots(4,1))) > int.inferior4<-c(mean(rep.inferior4[,1]), + mean(rep.inferior4[,1])-1.96*mean(na.omit(rep.inferior4[,2])), 3 + > > + + > > + + mean(rep.inferior4[,1])+1.96*mean(na.omit(rep.inferior4[,2]))) rep.medio4<-t(replicate(1000,boots(4,2))) int.medio4<-c(mean(rep.medio4[,1]), mean(rep.medio4[,1])-1.96*mean(na.omit(rep.medio4[,2])), mean(rep.medio4[,1])+1.96*mean(na.omit(rep.medio4[,2]))) rep.superior4<-t(replicate(1000,boots(4,3))) int.superior4<-c(mean(rep.superior4[,1]), mean(rep.superior4[,1])-1.96*mean(na.omit(rep.superior4[,2])), mean(rep.superior4[,1])+1.96*mean(na.omit(rep.superior4[,2]))) HC inf sup inferior 0.2912236 0.2373322 0.3451149 medio 0.3206657 0.2997713 0.3415601 superior 0.2967301 0.2445649 0.3488953 1.3 Lo que realmente interesaba a los autores del art´ıculo era estimar y estudiar la significaci´ on estad´ıstica de la diferencia de biodiversidades entre per´ıodos.
Calcula el intervalo de confianza bootstrap para las tres posibles diferencias de biodiversidades: 1. superior-medio, 2. superior-inferior y 3. medioinferior (no uses boot).
> boots2<-function(alpha,num1,num2){ + datos1<-sample(dinos.prob[,num1],replace=TRUE) + datos2<-sample(dinos.prob[,num2],replace=TRUE) + datos<-abs(datos2-datos1) + HC<-( 1/(alpha-1) ) * ( 1 - sum(datos^alpha)) + HC.se<-( 1/sqrt( length(datos) ) ) * abs( alpha/(alpha-1) ) * + sqrt( sum( datos^(2*alpha-1) ) - sum(datos^alpha)^2 ) + c(HC,HC.se) + } • α = 0.5 : > > + + > > + + > > + + rep.1.0.5<-t(replicate(1000,boots2(0.5,2,3))) int.1.0.5<-c(mean(rep.1.0.5[,1]), mean(rep.1.0.5[,1])-1.96*mean(na.omit(rep.1.0.5[,2])), mean(rep.1.0.5[,1])+1.96*mean(na.omit(rep.1.0.5[,2]))) rep.2.0.5<-t(replicate(1000,boots2(0.5,1,3))) int.2.0.5<-c(mean(rep.2.0.5[,1]), mean(rep.2.0.5[,1])-1.96*mean(na.omit(rep.2.0.5[,2])), mean(rep.2.0.5[,1])+1.96*mean(na.omit(rep.2.0.5[,2]))) rep.3.0.5<-t(replicate(1000,boots2(0.5,1,2))) int.3.0.5<-c(mean(rep.3.0.5[,1]), mean(rep.3.0.5[,1])-1.96*mean(na.omit(rep.3.0.5[,2])), mean(rep.3.0.5[,1])+1.96*mean(na.omit(rep.3.0.5[,2]))) HC dif inf sup superior-inferior 2.531134 1.469350 3.592919 superior-medio 2.705784 1.620444 3.791124 medio-inferior 2.644553 1.561208 3.727899 •α=2: 4 > > + + > > + + > > + + rep.1.2<-t(replicate(1000,boots2(2,2,3))) int.1.2<-c(mean(rep.1.2[,1]), mean(rep.1.2[,1])-1.96*mean(na.omit(rep.1.2[,2])), mean(rep.1.2[,1])+1.96*mean(na.omit(rep.1.2[,2]))) rep.2.2<-t(replicate(1000,boots2(2,1,3))) int.2.2<-c(mean(rep.2.2[,1]), mean(rep.2.2[,1])-1.96*mean(na.omit(rep.2.2[,2])), mean(rep.2.2[,1])+1.96*mean(na.omit(rep.2.2[,2]))) rep.3.2<-t(replicate(1000,boots2(2,1,2))) int.3.2<-c(mean(rep.3.2[,1]), mean(rep.3.2[,1])-1.96*mean(na.omit(rep.3.2[,2])), mean(rep.3.2[,1])+1.96*mean(na.omit(rep.3.2[,2]))) HC dif inf sup superior-inferior 0.6466744 0.4543716 0.8389771 superior-medio 0.4987598 0.2580260 0.7394936 medio-inferior 0.5912777 0.3826790 0.7998764 •α=4: > > + + > > + + > > + + rep.1.4<-t(replicate(1000,boots2(4,2,3))) int.1.4<-c(mean(rep.1.4[,1]), mean(rep.1.4[,1])-1.96*mean(na.omit(rep.1.4[,2])), mean(rep.1.4[,1])+1.96*mean(na.omit(rep.1.4[,2]))) rep.2.4<-t(replicate(1000,boots2(4,1,3))) int.2.4<-c(mean(rep.2.4[,1]), mean(rep.2.4[,1])-1.96*mean(na.omit(rep.2.4[,2])), mean(rep.2.4[,1])+1.96*mean(na.omit(rep.2.4[,2]))) rep.3.4<-t(replicate(1000,boots2(4,1,2))) int.3.4<-c(mean(rep.3.4[,1]), mean(rep.3.4[,1])-1.96*mean(na.omit(rep.3.4[,2])), mean(rep.3.4[,1])+1.96*mean(na.omit(rep.3.4[,2]))) HC dif inf sup superior-inferior 0.3131291 0.2808306 0.3454277 superior-medio 0.2944692 0.2459397 0.3429988 medio-inferior 0.3088197 0.2785401 0.3390993 1.4 Replica b. y c. usando la librer´ıa boot.
> require(boot) • Replica del apartado b. : • α = 0.5: > alpha<-0.5 > para.boot<-function(datos,id){ + ( 1/(alpha-1) ) * ( 1 - sum(datos[id]^alpha)) + } 5 > > + + > > + + > > + + boot.inferior<-boot(dinos.prob[,1],statistic=para.boot,R=1000) int.inferior0.5<-c(boot.inferior$t0, boot.inferior$t0-1.96*sd(boot.inferior$t), boot.inferior$t0+1.96*sd(boot.inferior$t)) boot.medio<-boot(dinos.prob[,2],statistic=para.boot,R=1000) int.medio0.5<-c(boot.medio$t0, boot.medio$t0-1.96*sd(boot.medio$t), boot.medio$t0+1.96*sd(boot.medio$t)) boot.superior<-boot(dinos.prob[,3],statistic=para.boot,R=1000) int.superior0.5<-c(boot.superior$t0, boot.superior$t0-1.96*sd(boot.superior$t), boot.superior$t0+1.96*sd(boot.superior$t)) HC inf sup inferior 1.831826 -0.5620279 4.225680 medio 2.359374 0.6437013 4.075047 superior 2.230405 0.2064264 4.254385 • α = 2: > alpha<-2 > para.boot<-function(datos,id){ + ( 1/(alpha-1) ) * ( 1 - sum(datos[id]^alpha)) + } > > + + > > + + > > + + boot.inferior<-boot(dinos.prob[,1],statistic=para.boot,R=1000) int.inferior2<-c(boot.inferior$t0, boot.inferior$t0-1.96*sd(boot.inferior$t), boot.inferior$t0+1.96*sd(boot.inferior$t)) boot.medio<-boot(dinos.prob[,2],statistic=para.boot,R=1000) int.medio2<-c(boot.medio$t0, boot.medio$t0-1.96*sd(boot.medio$t), boot.medio$t0+1.96*sd(boot.medio$t)) boot.superior<-boot(dinos.prob[,3],statistic=para.boot,R=1000) int.superior2<-c(boot.superior$t0, boot.superior$t0-1.96*sd(boot.superior$t), boot.superior$t0+1.96*sd(boot.superior$t)) HC inf sup inferior 0.5556877 -0.05250204 1.1638774 medio 0.7000000 0.40828545 0.9917146 superior 0.6097337 0.02252757 1.1969398 • α = 4: > alpha<-4 > para.boot<-function(datos,id){ + ( 1/(alpha-1) ) * ( 1 - sum(datos[id]^alpha)) + } > boot.inferior<-boot(dinos.prob[,1],statistic=para.boot,R=1000) > int.inferior4<-c(boot.inferior$t0, 6 + + > > + + > > + + > boot.inferior$t0-1.96*sd(boot.inferior$t), boot.inferior$t0+1.96*sd(boot.inferior$t)) boot.medio<-boot(dinos.prob[,2],statistic=para.boot,R=1000) int.medio4<-c(boot.medio$t0, boot.medio$t0-1.96*sd(boot.medio$t), boot.medio$t0+1.96*sd(boot.medio$t)) boot.superior<-boot(dinos.prob[,3],statistic=para.boot,R=1000) int.superior4<-c(boot.superior$t0, boot.superior$t0-1.96*sd(boot.superior$t), boot.superior$t0+1.96*sd(boot.superior$t)) HC inf sup inferior 0.2928517 0.2290267 0.3566766 medio 0.3204989 0.3057345 0.3352632 superior 0.2960006 0.2326408 0.3593604 • Replica del apartado b. : • α = 0.5: > alpha<-0.5 > para.boot<-function(datos,id){ + ( 1/(alpha-1) ) * ( 1 - sum(datos[id]^alpha)) + } > + > + + > + > + + > + > + + boot.1<-boot(abs(dinos.prob[,3]-dinos.prob[,2]), statistic=para.boot,R=1000) int.1.0.5<-c(boot.1$t0, boot.1$t0-1.96*sd(boot.1$t), boot.1$t0+1.96*sd(boot.1$t)) boot.2<-boot(abs(dinos.prob[,3]-dinos.prob[,1]), statistic=para.boot,R=1000) int.2.0.5<-c(boot.2$t0, boot.2$t0-1.96*sd(boot.2$t), boot.2$t0+1.96*sd(boot.2$t)) boot.3<-boot(abs(dinos.prob[,2]-dinos.prob[,1]), statistic=para.boot,R=1000) int.3.0.5<-c(boot.3$t0, boot.3$t0-1.96*sd(boot.3$t), boot.3$t0+1.96*sd(boot.3$t)) HC inf sup superior-inferior 1.4795588 0.4993838 2.459734 superior-medio 0.1892394 -0.8691311 1.247610 medio-inferior 0.6739863 -0.5251295 1.873102 • α = 2: > alpha<-2 > para.boot<-function(datos,id){ 7 + ( 1/(alpha-1) ) * ( 1 - sum(datos[id]^alpha)) + } > + > + + > + > + + > + > + + boot.1<-boot(abs(dinos.prob[,3]-dinos.prob[,2]), statistic=para.boot,R=1000) int.1.2<-c(boot.1$t0, boot.1$t0-1.96*sd(boot.1$t), boot.1$t0+1.96*sd(boot.1$t)) boot.2<-boot(abs(dinos.prob[,3]-dinos.prob[,1]), statistic=para.boot,R=1000) int.2.2<-c(boot.2$t0, boot.2$t0-1.96*sd(boot.2$t), boot.2$t0+1.96*sd(boot.2$t)) boot.3<-boot(abs(dinos.prob[,2]-dinos.prob[,1]), statistic=para.boot,R=1000) int.3.2<-c(boot.3$t0, boot.3$t0-1.96*sd(boot.3$t), boot.3$t0+1.96*sd(boot.3$t)) HC inf sup superior-inferior 0.9222878 0.8525348 0.9920409 superior-medio 0.9762887 0.9502515 1.0023258 medio-inferior 0.9466565 0.8836802 1.0096328 • α = 4: > alpha<-4 > para.boot<-function(datos,id){ + ( 1/(alpha-1) ) * ( 1 - sum(datos[id]^alpha)) + } > + > + + > + > + + > + > + + boot.1<-boot(abs(dinos.prob[,3]-dinos.prob[,2]), statistic=para.boot,R=1000) int.1.4<-c(boot.1$t0, boot.1$t0-1.96*sd(boot.1$t), boot.1$t0+1.96*sd(boot.1$t)) boot.2<-boot(abs(dinos.prob[,3]-dinos.prob[,1]), statistic=para.boot,R=1000) int.2.4<-c(boot.2$t0, boot.2$t0-1.96*sd(boot.2$t), boot.2$t0+1.96*sd(boot.2$t)) boot.3<-boot(abs(dinos.prob[,2]-dinos.prob[,1]), statistic=para.boot,R=1000) int.3.4<-c(boot.3$t0, boot.3$t0-1.96*sd(boot.3$t), boot.3$t0+1.96*sd(boot.3$t)) HC inf sup superior-inferior 0.3326068 0.3316925 0.3335211 superior-medio 0.3332462 0.3331288 0.3333635 medio-inferior 0.3327855 0.3319022 0.3336688 8 1.5 Bas´ andonos en los resultados de bootstrap, ¿qu´ e crees que pas´ o con los dinosaurios? En todos los casos se ve que los tres per´ıodos dan estimaciones y errores muy parecidos. De hecho, en los estudios sobre las diferencias entre per´ıodos, queda claro que las diferencias entre ellos son iguales, indistintamente de cu´ ales se comparen.
En conclusi´ on, se podr´ıa concluir que la desaparici´on de los dinosaurios fue gradual, ya que en todos los per´ıodos desaparecieron a la misma velocidad.
9 ...

Tags: