Tarea uno
Rocío Mendieta Oviedo
22/8/2022
Tarea uno de la materia de REGRESIÓN APLICADA A LAS CIENCIAS AMBIENTALES impartida por el Dr. Jorge Méndez González.
Variables generadas
a) Con los datos data(trees) del paquete de r se generaron las siguientes variables: Girth2, Height2, GirthxHeight, Girth2*Height,con la finalidad de tener un total de 8 variables. Para realizar las trasformaciones se apoyó del script https://rpubs.com/jorge_mendez/601482 en la sección donde dice “Hacer transformaciones y renombrar variables”.
setwd("C:/Qto_semestre/Regresión") # Ajustar nuestro directorio
data ("trees") # Conjunto de datos
Datos<-trees
head(Datos) # Primeros datos## Girth Height Volume
## 1 8.3 70 10.3
## 2 8.6 65 10.3
## 3 8.8 63 10.2
## 4 10.5 72 16.4
## 5 10.7 81 18.8
## 6 10.8 83 19.7str(Datos) # ver la estructura de los datos## 'data.frame': 31 obs. of 3 variables:
## $ Girth : num 8.3 8.6 8.8 10.5 10.7 10.8 11 11 11.1 11.2 ...
## $ Height: num 70 65 63 72 81 83 66 75 80 75 ...
## $ Volume: num 10.3 10.3 10.2 16.4 18.8 19.7 15.6 18.2 22.6 19.9 ...attach(Datos) # Comando crucial
Girth2 <- Girth * Girth
Height2 <- Height* Height
GirthxHeight <- Girth*Height
Girth2xHeight<-Girth2*Height
Datos<-data.frame(Girth, Height, Volume, Girth2, Height2, GirthxHeight, Girth2xHeight )
attach(Datos) #La base de datos prelimina## The following objects are masked _by_ .GlobalEnv:
##
## Girth2, Girth2xHeight, GirthxHeight, Height2## The following objects are masked from Datos (pos = 3):
##
## Girth, Height, Volumeb) Después, de la tabla obtenida con las ocho variables se obtuvo una muestra aleatoria de 25 datos, para hacer los análisis estadísticos de correlación de pearson. Antes de ello se istalaron las library necesarias para los analisis estadisticos
library(sampler)
library(tidyverse)## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
## ✔ ggplot2 3.3.6 ✔ purrr 0.3.4
## ✔ tibble 3.1.8 ✔ dplyr 1.0.9
## ✔ tidyr 1.2.0 ✔ stringr 1.4.1
## ✔ readr 2.1.2 ✔ forcats 0.5.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()library(dplyr)
library(ggplot2)
library(GGally)## Registered S3 method overwritten by 'GGally':
## method from
## +.gg ggplot2library(Hmisc)## Loading required package: lattice
## Loading required package: survival
## Loading required package: Formula
##
## Attaching package: 'Hmisc'
##
## The following objects are masked from 'package:dplyr':
##
## src, summarize
##
## The following objects are masked from 'package:base':
##
## format.pval, unitslibrary(corrplot)## corrplot 0.92 loadedlibrary(DescTools)##
## Attaching package: 'DescTools'
##
## The following objects are masked from 'package:Hmisc':
##
## %nin%, Label, Mean, Quantilelibrary(tidyr)
library(pastecs)##
## Attaching package: 'pastecs'
##
## The following objects are masked from 'package:dplyr':
##
## first, last
##
## The following object is masked from 'package:tidyr':
##
## extractlibrary(nortest)
# library(normtest)
library(mvnormtest)Muestra aleatoria de 25 datos
obtenemos la muestra y exportamos el archivo en csv
Datos_m<-rsamp(Datos, n=25, rep=FALSE); Datos_m # Obtenemos la muestra## Girth Height Volume Girth2 Height2 GirthxHeight Girth2xHeight
## 1 17.9 80 58.3 320.41 6400 1432.0 25632.80
## 2 12.0 75 19.1 144.00 5625 900.0 10800.00
## 3 8.3 70 10.3 68.89 4900 581.0 4822.30
## 4 14.5 74 36.3 210.25 5476 1073.0 15558.50
## 5 18.0 80 51.5 324.00 6400 1440.0 25920.00
## 6 8.6 65 10.3 73.96 4225 559.0 4807.40
## 7 11.1 80 22.6 123.21 6400 888.0 9856.80
## 8 11.3 79 24.2 127.69 6241 892.7 10087.51
## 9 11.2 75 19.9 125.44 5625 840.0 9408.00
## 10 14.2 80 31.7 201.64 6400 1136.0 16131.20
## 11 13.7 71 25.7 187.69 5041 972.7 13325.99
## 12 12.9 85 33.8 166.41 7225 1096.5 14144.85
## 13 14.0 78 34.5 196.00 6084 1092.0 15288.00
## 14 8.8 63 10.2 77.44 3969 554.4 4878.72
## 15 10.5 72 16.4 110.25 5184 756.0 7938.00
## 16 10.8 83 19.7 116.64 6889 896.4 9681.12
## 17 20.6 87 77.0 424.36 7569 1792.2 36919.32
## 18 13.8 64 24.9 190.44 4096 883.2 12188.16
## 19 18.0 80 51.0 324.00 6400 1440.0 25920.00
## 20 16.0 72 38.3 256.00 5184 1152.0 18432.00
## 21 11.0 66 15.6 121.00 4356 726.0 7986.00
## 22 16.3 77 42.6 265.69 5929 1255.1 20458.13
## 23 13.3 86 27.4 176.89 7396 1143.8 15212.54
## 24 11.7 69 21.3 136.89 4761 807.3 9445.41
## 25 11.4 76 21.4 129.96 5776 866.4 9876.96# write.csv(Datos_m, file = "Muestra.csv") # Exportamos el archivo.Correlación
c) Se aplico correlación de Pearson con ayuda de este video (https://www.youtube.com/watch?v=qFIXTihSZ2w) para analizar cuál de todas las variables posee mayor correlación con la variable volumen.
dat=read.csv("C:/Qto_semestre/Regresión/Muestra.csv")
attach(dat);dat## The following objects are masked _by_ .GlobalEnv:
##
## Girth2, GirthxHeight, Height2## The following objects are masked from Datos (pos = 23):
##
## Girth, Girth2, GirthxHeight, Height, Height2## The following objects are masked from Datos (pos = 24):
##
## Girth, Height## X Girth Height Volumennn Girth2 Height2 GirthxHeight diam2xalt
## 1 1 17.5 82 55.7 306.25 6724 1435.0 25112.50
## 2 2 12.0 75 19.1 144.00 5625 900.0 10800.00
## 3 3 16.0 72 38.3 256.00 5184 1152.0 18432.00
## 4 4 11.1 80 22.6 123.21 6400 888.0 9856.80
## 5 5 13.7 71 25.7 187.69 5041 972.7 13325.99
## 6 6 10.8 83 19.7 116.64 6889 896.4 9681.12
## 7 7 11.3 79 24.2 127.69 6241 892.7 10087.51
## 8 8 11.4 76 21.0 129.96 5776 866.4 9876.96
## 9 9 11.4 76 21.4 129.96 5776 866.4 9876.96
## 10 10 18.0 80 51.0 324.00 6400 1440.0 25920.00
## 11 11 13.3 86 27.4 176.89 7396 1143.8 15212.54
## 12 12 11.7 69 21.3 136.89 4761 807.3 9445.41
## 13 13 12.9 74 22.2 166.41 5476 954.6 12314.34
## 14 14 11.0 66 15.6 121.00 4356 726.0 7986.00
## 15 15 17.3 81 55.4 299.29 6561 1401.3 24242.49
## 16 16 18.0 80 51.5 324.00 6400 1440.0 25920.00
## 17 17 13.8 64 24.9 190.44 4096 883.2 12188.16
## 18 18 20.6 87 77.0 424.36 7569 1792.2 36919.32
## 19 19 14.0 78 34.5 196.00 6084 1092.0 15288.00
## 20 20 16.3 77 42.6 265.69 5929 1255.1 20458.13
## 21 21 12.9 85 33.8 166.41 7225 1096.5 14144.85
## 22 22 8.3 70 10.3 68.89 4900 581.0 4822.30
## 23 23 17.9 80 58.3 320.41 6400 1432.0 25632.80
## 24 24 14.2 80 31.7 201.64 6400 1136.0 16131.20
## 25 25 10.7 81 18.8 114.49 6561 866.7 9273.69####################################### Analisis exploratorio de los datos en plots
library(PerformanceAnalytics)## Loading required package: xts## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric##
## Attaching package: 'xts'## The following objects are masked from 'package:pastecs':
##
## first, last## The following objects are masked from 'package:dplyr':
##
## first, last##
## Attaching package: 'PerformanceAnalytics'## The following object is masked from 'package:graphics':
##
## legend# Plot correlacion todas las variables -------------------------------------------
chart.Correlation(dat[,2:8], # elegir las variables
method="pearson", # Puede usar otro tipo de correlacion
histogram=TRUE,
pch=16)## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
library(correlation)
correlation(dat, method = "pearson")## # Correlation Matrix (pearson-method)
##
## Parameter1 | Parameter2 | r | 95% CI | t(23) | p
## -----------------------------------------------------------------------
## X | Girth | 0.09 | [-0.32, 0.47] | 0.42 | > .999
## X | Height | 0.09 | [-0.32, 0.47] | 0.41 | > .999
## X | Volumennn | 0.14 | [-0.27, 0.51] | 0.69 | > .999
## X | Girth2 | 0.11 | [-0.30, 0.48] | 0.52 | > .999
## X | Height2 | 0.09 | [-0.31, 0.47] | 0.46 | > .999
## X | GirthxHeight | 0.11 | [-0.30, 0.49] | 0.54 | > .999
## X | diam2xalt | 0.12 | [-0.28, 0.50] | 0.60 | > .999
## Girth | Height | 0.39 | [-0.01, 0.68] | 2.00 | 0.501
## Girth | Volumennn | 0.96 | [ 0.91, 0.98] | 16.84 | < .001***
## Girth | Girth2 | 0.99 | [ 0.99, 1.00] | 43.84 | < .001***
## Girth | Height2 | 0.39 | [-0.01, 0.68] | 2.03 | 0.501
## Girth | GirthxHeight | 0.97 | [ 0.92, 0.98] | 17.65 | < .001***
## Girth | diam2xalt | 0.98 | [ 0.95, 0.99] | 22.38 | < .001***
## Height | Volumennn | 0.52 | [ 0.16, 0.76] | 2.92 | 0.107
## Height | Girth2 | 0.40 | [ 0.00, 0.69] | 2.08 | 0.501
## Height | Height2 | 1.00 | [ 1.00, 1.00] | 107.10 | < .001***
## Height | GirthxHeight | 0.61 | [ 0.28, 0.81] | 3.67 | 0.020*
## Height | diam2xalt | 0.51 | [ 0.14, 0.75] | 2.84 | 0.113
## Volumennn | Girth2 | 0.98 | [ 0.94, 0.99] | 21.11 | < .001***
## Volumennn | Height2 | 0.53 | [ 0.16, 0.76] | 2.96 | 0.104
## Volumennn | GirthxHeight | 0.98 | [ 0.95, 0.99] | 22.49 | < .001***
## Volumennn | diam2xalt | 0.99 | [ 0.97, 0.99] | 29.92 | < .001***
## Girth2 | Height2 | 0.40 | [ 0.01, 0.69] | 2.11 | 0.501
## Girth2 | GirthxHeight | 0.97 | [ 0.93, 0.99] | 18.69 | < .001***
## Girth2 | diam2xalt | 0.99 | [ 0.98, 1.00] | 34.24 | < .001***
## Height2 | GirthxHeight | 0.61 | [ 0.29, 0.81] | 3.71 | 0.019*
## Height2 | diam2xalt | 0.51 | [ 0.15, 0.76] | 2.88 | 0.110
## GirthxHeight | diam2xalt | 0.99 | [ 0.97, 1.00] | 31.37 | < .001***
##
## p-value adjustment method: Holm (1979)
## Observations: 25Analizando los resultados del inciso c) la variable que tiene mayor correlación con volume (Volumennn) es Girth2*Height (diam2xalt), ya que se tuvo la variable con mayor correlación se analizó la normalidad entre ambas variables para aplicar correlación de pearson ya que si la muestra no proviene de una distribución normal no se puede realizar dicha correlación.
Normalidad
d) Se Verifico que se cumpliera el supuesto de normalidad en ambas variables, Volumennn y diam2xalt
Hipótesis
H0: La muestra proviene de una distribución normal.
H1: La muestra no proviene de una distribución normal.
Nivel de Significancia
El nivel de significancia que se trabajará es de 0.05. Alfa=0.05
Criterio de Decisión
Si P < Alfa Se rechaza Ho
Si p >= Alfa No se rechaza Ho
############ Volumennn
ad.test(Volumennn) # Prueba de Anderson-Darling ##
## Anderson-Darling normality test
##
## data: Volumennn
## A = 1.1151, p-value = 0.005186cvm.test(Volumennn ) # Prueba de Cramer-von Mises, muestras pequeñas##
## Cramer-von Mises normality test
##
## data: Volumennn
## W = 0.20188, p-value = 0.004535lillie.test(Volumennn ) # Prueba de Lilliefors (Kolmogorov-Smirnov) muestras grandes##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: Volumennn
## D = 0.19089, p-value = 0.01936sf.test(Volumennn ) # Prueba de Shapiro-Francia muestras pequeñas##
## Shapiro-Francia normality test
##
## data: Volumennn
## W = 0.89068, p-value = 0.01378shapiro.test(Volumennn ) # Prueba de normalidad##
## Shapiro-Wilk normality test
##
## data: Volumennn
## W = 0.89104, p-value = 0.01177############ diam2xalt
ad.test(diam2xalt) # Prueba de Anderson-Darling ##
## Anderson-Darling normality test
##
## data: diam2xalt
## A = 1.0963, p-value = 0.005789cvm.test(diam2xalt) # Prueba de Cramer-von Mises, muestras pequeñas##
## Cramer-von Mises normality test
##
## data: diam2xalt
## W = 0.19037, p-value = 0.006382lillie.test(diam2xalt) # Prueba de Lilliefors (Kolmogorov-Smirnov) muestras grandes##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: diam2xalt
## D = 0.16218, p-value = 0.08865sf.test(diam2xalt) # Prueba de Shapiro-Francia muestras pequeñas##
## Shapiro-Francia normality test
##
## data: diam2xalt
## W = 0.88599, p-value = 0.01132shapiro.test(diam2xalt) # Prueba de n##
## Shapiro-Wilk normality test
##
## data: diam2xalt
## W = 0.88865, p-value = 0.01047Visiblemente los valores de p-value son menores a 0.05 en ambas variables aplicando diferentes pruebas de normalidad por lo tanto se rechaza hipótesis nula, entonces es necesaria una trasformación, se opto por Box Cox para que las variables cumplan el supuesto de normalidad.
Trasformacióm Box Cox
e) Se realizó transformación box Cox.
# Transformación Box - Cox _Volumennn
library(car)## Loading required package: carData##
## Attaching package: 'car'## The following object is masked from 'package:DescTools':
##
## Recode## The following object is masked from 'package:dplyr':
##
## recode## The following object is masked from 'package:purrr':
##
## somepowerTransform(Volumennn )## Estimated transformation parameter
## Volumennn
## -0.115556y_1=Volumennn ^ -0.115556
hist(y_1) 
shapiro.test(y_1)##
## Shapiro-Wilk normality test
##
## data: y_1
## W = 0.96487, p-value = 0.5197ad.test(y_1) # Prueba de Anderson-Darling ##
## Anderson-Darling normality test
##
## data: y_1
## A = 0.42461, p-value = 0.2934cvm.test(y_1) # Prueba de Cramer-von Mises, muestras pequeñas##
## Cramer-von Mises normality test
##
## data: y_1
## W = 0.069981, p-value = 0.2683lillie.test(y_1) # Prueba de Lilliefors (Kolmogorov-Smirnov) muestras grandes##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: y_1
## D = 0.11931, p-value = 0.475sf.test(y_1) # Prueba de Shapiro-Francia muestras pequeñas##
## Shapiro-Francia normality test
##
## data: y_1
## W = 0.96444, p-value = 0.4334shapiro.test(y_1) # Prueba de n##
## Shapiro-Wilk normality test
##
## data: y_1
## W = 0.96487, p-value = 0.5197# Transformación Box - Cox _diam2xalt
library(car)
powerTransform(diam2xalt)## Estimated transformation parameter
## diam2xalt
## -0.06135679x_1=diam2xalt^ -0.06135679
hist(x_1)
shapiro.test(x_1)##
## Shapiro-Wilk normality test
##
## data: x_1
## W = 0.95791, p-value = 0.3744ad.test(x_1) # Prueba de Anderson-Darling ##
## Anderson-Darling normality test
##
## data: x_1
## A = 0.48972, p-value = 0.2015cvm.test(x_1) # Prueba de Cramer-von Mises, muestras pequeñas##
## Cramer-von Mises normality test
##
## data: x_1
## W = 0.073366, p-value = 0.2419lillie.test(x_1) # Prueba de Lilliefors (Kolmogorov-Smirnov) muestras grandes##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: x_1
## D = 0.11646, p-value = 0.5145sf.test(x_1) # Prueba de Shapiro-Francia muestras pequeñas##
## Shapiro-Francia normality test
##
## data: x_1
## W = 0.95535, p-value = 0.2808shapiro.test(x_1) # Prueba de n##
## Shapiro-Wilk normality test
##
## data: x_1
## W = 0.95791, p-value = 0.3744ya que se realizó dicha transformación se obtuvieron valores de p-value mayores a 0.05, entonces no se rechaza hipótesis nula y ahora si se puede realizar correlación de pearson considerando que ambas variables ahora si provienen de una muestra con distribución normal.
Correlación de pearson
f) Una vez que se cumplio la normalidad se realizó la prueba de significancia estadística (95%) de de la correlación de pearson.
data<-data.frame( y_1, x_1)
attach(data) ## The following objects are masked _by_ .GlobalEnv:
##
## x_1, y_1Pearson <- cor_test(data, "y_1", "x_1", method = "pearson"); Pearson # Con normalidad## Parameter1 | Parameter2 | r | 95% CI | t(23) | p
## -----------------------------------------------------------------
## y_1 | x_1 | 0.98 | [0.96, 0.99] | 26.97 | < .001***
##
## Observations: 25Por lo tanto con la trasformación de Box Cox se tiene una probabilidad del 95% que las variables Volumennn y diam2xalt se distribuye normal con un p-value < 2.2e-16 y un valor de r de 0.9845582. cumpliendo la significancia estadística la correlación entre las variables antes descritas se procedió a crear los intervalos de confianza (IC) y un Bootstrap cabe mencionar que si no se cumplía la significancia estadística no serian correctos los IC y el Bootstrap
Intervalos de confianza de la correlacion
g) Se obtuvieron los intervalos de confianza de la correlación de Pearson. ####################################### Intervalos de confianza de la correlacion
Pearson <- cor_test(data, "y_1", "x_1", method = "pearson"); Pearson## Parameter1 | Parameter2 | r | 95% CI | t(23) | p
## -----------------------------------------------------------------
## y_1 | x_1 | 0.98 | [0.96, 0.99] | 26.97 | < .001***
##
## Observations: 25Se tiene un intervalo de confianza de [0.96, 0.99] para la correlación de dichas variables.
Bootstrap correlacion
h) finalmente se Realizó un bootstrap de la correlación de Pearson con 500 repeticiones.
Bootstrap correlacion
library(boot)##
## Attaching package: 'boot'## The following object is masked from 'package:car':
##
## logit## The following object is masked from 'package:survival':
##
## aml## The following object is masked from 'package:lattice':
##
## melanomadata<-data.frame(y_1, x_1 )
Mi_corr<-cor.test(y_1, x_1) # Ya esta nornalizada, por eso usamos Pearson
str(Mi_corr)## List of 9
## $ statistic : Named num 27
## ..- attr(*, "names")= chr "t"
## $ parameter : Named int 23
## ..- attr(*, "names")= chr "df"
## $ p.value : num 6.54e-19
## $ estimate : Named num 0.985
## ..- attr(*, "names")= chr "cor"
## $ null.value : Named num 0
## ..- attr(*, "names")= chr "correlation"
## $ alternative: chr "two.sided"
## $ method : chr "Pearson's product-moment correlation"
## $ data.name : chr "y_1 and x_1"
## $ conf.int : num [1:2] 0.965 0.993
## ..- attr(*, "conf.level")= num 0.95
## - attr(*, "class")= chr "htest"Mi_corr$ estimate #Aqui vemos el valor de r## cor
## 0.9845582str(data)## 'data.frame': 25 obs. of 2 variables:
## $ y_1: num 0.628 0.711 0.656 0.697 0.687 ...
## $ x_1: num 0.537 0.566 0.547 0.569 0.558 ...set.seed(1)
b3 <- boot(data,
statistic = function(data, i) {
cor(data[i, "y_1"], data[i, "x_1"], method='pearson')
},
R = 500
)
b3##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = data, statistic = function(data, i) {
## cor(data[i, "y_1"], data[i, "x_1"], method = "pearson")
## }, R = 500)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 0.9845582 -0.001337225 0.007712174boot.ci(b3, type = c("norm", "basic", "perc", "bca"))## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 500 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = b3, type = c("norm", "basic", "perc", "bca"))
##
## Intervals :
## Level Normal Basic
## 95% ( 0.9708, 1.0010 ) ( 0.9759, 1.0035 )
##
## Level Percentile BCa
## 95% ( 0.9656, 0.9932 ) ( 0.9656, 0.9932 )
## Calculations and Intervals on Original Scaleplot(density(b3$t))
abline(v = Mi_corr$ estimate, lty = "dashed", col = "red") # Traemos el velor desde la corr
b3$t # veamos todos los resutados del bootstrapped## [,1]
## [1,] 0.9864679
## [2,] 0.9792835
## [3,] 0.9907636
## [4,] 0.9865235
## [5,] 0.9860830
## [6,] 0.9869526
## [7,] 0.9900609
## [8,] 0.9883850
## [9,] 0.9873034
## [10,] 0.9806165
## [11,] 0.9704991
## [12,] 0.9794573
## [13,] 0.9815366
## [14,] 0.9820594
## [15,] 0.9859439
## [16,] 0.9834093
## [17,] 0.9882646
## [18,] 0.9784452
## [19,] 0.9800394
## [20,] 0.9767491
## [21,] 0.9890766
## [22,] 0.9800252
## [23,] 0.9685309
## [24,] 0.9879794
## [25,] 0.9765526
## [26,] 0.9842044
## [27,] 0.9860087
## [28,] 0.9824461
## [29,] 0.9809522
## [30,] 0.9863616
## [31,] 0.9764422
## [32,] 0.9774163
## [33,] 0.9817254
## [34,] 0.9752856
## [35,] 0.9728600
## [36,] 0.9870425
## [37,] 0.9837736
## [38,] 0.9712038
## [39,] 0.9911194
## [40,] 0.9774832
## [41,] 0.9880089
## [42,] 0.9837906
## [43,] 0.9912561
## [44,] 0.9889168
## [45,] 0.9907847
## [46,] 0.9847716
## [47,] 0.9909470
## [48,] 0.9789277
## [49,] 0.9900387
## [50,] 0.9770805
## [51,] 0.9858427
## [52,] 0.9821416
## [53,] 0.9900574
## [54,] 0.9690063
## [55,] 0.9799891
## [56,] 0.9774956
## [57,] 0.9821977
## [58,] 0.9875638
## [59,] 0.9916158
## [60,] 0.9820549
## [61,] 0.9701176
## [62,] 0.9848207
## [63,] 0.9925879
## [64,] 0.9850320
## [65,] 0.9839902
## [66,] 0.9889212
## [67,] 0.9873977
## [68,] 0.9927222
## [69,] 0.9815068
## [70,] 0.9797303
## [71,] 0.9842961
## [72,] 0.9867697
## [73,] 0.9870261
## [74,] 0.9841977
## [75,] 0.9831734
## [76,] 0.9853016
## [77,] 0.9862689
## [78,] 0.9781457
## [79,] 0.9836668
## [80,] 0.9871505
## [81,] 0.9898694
## [82,] 0.9801199
## [83,] 0.9835249
## [84,] 0.9864735
## [85,] 0.9848388
## [86,] 0.9869067
## [87,] 0.9836841
## [88,] 0.9869525
## [89,] 0.9836939
## [90,] 0.9904228
## [91,] 0.9895791
## [92,] 0.9890103
## [93,] 0.9803200
## [94,] 0.9609764
## [95,] 0.9777815
## [96,] 0.9737311
## [97,] 0.9765066
## [98,] 0.9844651
## [99,] 0.9891734
## [100,] 0.9896863
## [101,] 0.9887903
## [102,] 0.9820095
## [103,] 0.9882174
## [104,] 0.9691244
## [105,] 0.9835564
## [106,] 0.9595129
## [107,] 0.9796085
## [108,] 0.9848177
## [109,] 0.9819562
## [110,] 0.9863156
## [111,] 0.9881551
## [112,] 0.9781720
## [113,] 0.9458501
## [114,] 0.9839783
## [115,] 0.9885285
## [116,] 0.9831711
## [117,] 0.9900066
## [118,] 0.9679045
## [119,] 0.9789433
## [120,] 0.9865355
## [121,] 0.9888432
## [122,] 0.9885336
## [123,] 0.9801784
## [124,] 0.9885850
## [125,] 0.9821309
## [126,] 0.9664905
## [127,] 0.9828882
## [128,] 0.9897236
## [129,] 0.9927014
## [130,] 0.9922527
## [131,] 0.9766920
## [132,] 0.9910212
## [133,] 0.9877460
## [134,] 0.9739331
## [135,] 0.9851099
## [136,] 0.9884294
## [137,] 0.9833357
## [138,] 0.9761388
## [139,] 0.9888527
## [140,] 0.9877003
## [141,] 0.9698340
## [142,] 0.9877919
## [143,] 0.9886508
## [144,] 0.9720532
## [145,] 0.9840361
## [146,] 0.9920791
## [147,] 0.9905997
## [148,] 0.9873037
## [149,] 0.9796185
## [150,] 0.9910070
## [151,] 0.9872218
## [152,] 0.9898775
## [153,] 0.9845621
## [154,] 0.9926549
## [155,] 0.9853411
## [156,] 0.9958222
## [157,] 0.9888224
## [158,] 0.9794728
## [159,] 0.9914676
## [160,] 0.9734417
## [161,] 0.9874399
## [162,] 0.9667698
## [163,] 0.9804119
## [164,] 0.9921075
## [165,] 0.9872708
## [166,] 0.9930574
## [167,] 0.9818036
## [168,] 0.9934375
## [169,] 0.9869983
## [170,] 0.9896225
## [171,] 0.9895949
## [172,] 0.9870567
## [173,] 0.9797336
## [174,] 0.9851938
## [175,] 0.9829246
## [176,] 0.9936555
## [177,] 0.9877974
## [178,] 0.9815504
## [179,] 0.9862241
## [180,] 0.9822731
## [181,] 0.9340007
## [182,] 0.9860206
## [183,] 0.9871142
## [184,] 0.9827809
## [185,] 0.9813356
## [186,] 0.9878227
## [187,] 0.9911538
## [188,] 0.9825994
## [189,] 0.9874054
## [190,] 0.9945549
## [191,] 0.9811268
## [192,] 0.9808472
## [193,] 0.9710154
## [194,] 0.9849265
## [195,] 0.9890356
## [196,] 0.9800607
## [197,] 0.9674740
## [198,] 0.9776713
## [199,] 0.9852945
## [200,] 0.9917349
## [201,] 0.9860265
## [202,] 0.9740164
## [203,] 0.9778946
## [204,] 0.9920231
## [205,] 0.9856767
## [206,] 0.9895826
## [207,] 0.9739121
## [208,] 0.9880494
## [209,] 0.9879055
## [210,] 0.9841412
## [211,] 0.9872867
## [212,] 0.9870035
## [213,] 0.9869478
## [214,] 0.9866831
## [215,] 0.9835514
## [216,] 0.9889829
## [217,] 0.9920477
## [218,] 0.9765110
## [219,] 0.9773773
## [220,] 0.9837025
## [221,] 0.9805446
## [222,] 0.9890367
## [223,] 0.9898242
## [224,] 0.9907184
## [225,] 0.9834739
## [226,] 0.9873997
## [227,] 0.9747536
## [228,] 0.9806587
## [229,] 0.9892174
## [230,] 0.9885145
## [231,] 0.9874287
## [232,] 0.9740492
## [233,] 0.9813051
## [234,] 0.9916197
## [235,] 0.9897496
## [236,] 0.9890943
## [237,] 0.9924195
## [238,] 0.9857368
## [239,] 0.9813315
## [240,] 0.9545879
## [241,] 0.9749617
## [242,] 0.9829654
## [243,] 0.9833155
## [244,] 0.9826140
## [245,] 0.9824183
## [246,] 0.9856508
## [247,] 0.9901048
## [248,] 0.9721886
## [249,] 0.9933398
## [250,] 0.9761008
## [251,] 0.9838840
## [252,] 0.9748777
## [253,] 0.9898008
## [254,] 0.9785190
## [255,] 0.9658350
## [256,] 0.9814501
## [257,] 0.9895409
## [258,] 0.9876556
## [259,] 0.9805612
## [260,] 0.9884018
## [261,] 0.9842507
## [262,] 0.9917637
## [263,] 0.9899232
## [264,] 0.9852329
## [265,] 0.9814153
## [266,] 0.9830876
## [267,] 0.9753711
## [268,] 0.9827044
## [269,] 0.9862291
## [270,] 0.9807736
## [271,] 0.9919837
## [272,] 0.9809864
## [273,] 0.9819731
## [274,] 0.9785403
## [275,] 0.9803478
## [276,] 0.9751081
## [277,] 0.9744414
## [278,] 0.9814345
## [279,] 0.9669492
## [280,] 0.9882198
## [281,] 0.9904739
## [282,] 0.9851520
## [283,] 0.9869330
## [284,] 0.9673865
## [285,] 0.9924278
## [286,] 0.9909705
## [287,] 0.9918117
## [288,] 0.9701022
## [289,] 0.9555095
## [290,] 0.9875529
## [291,] 0.9918963
## [292,] 0.9935940
## [293,] 0.9843682
## [294,] 0.9918958
## [295,] 0.9850091
## [296,] 0.9889610
## [297,] 0.9765109
## [298,] 0.9714010
## [299,] 0.9454864
## [300,] 0.9883810
## [301,] 0.9724000
## [302,] 0.9845381
## [303,] 0.9900613
## [304,] 0.9878991
## [305,] 0.9846656
## [306,] 0.9862736
## [307,] 0.9788980
## [308,] 0.9716747
## [309,] 0.9862187
## [310,] 0.9856233
## [311,] 0.9831322
## [312,] 0.9787752
## [313,] 0.9838911
## [314,] 0.9774275
## [315,] 0.9917562
## [316,] 0.9868697
## [317,] 0.9793027
## [318,] 0.9717062
## [319,] 0.9833002
## [320,] 0.9916451
## [321,] 0.9710011
## [322,] 0.9756036
## [323,] 0.9853875
## [324,] 0.9832748
## [325,] 0.9904835
## [326,] 0.9819138
## [327,] 0.9711679
## [328,] 0.9829242
## [329,] 0.9818936
## [330,] 0.9901096
## [331,] 0.9895554
## [332,] 0.9790281
## [333,] 0.9873091
## [334,] 0.9931167
## [335,] 0.9887397
## [336,] 0.9736809
## [337,] 0.9921230
## [338,] 0.9649098
## [339,] 0.9890901
## [340,] 0.9923275
## [341,] 0.9899868
## [342,] 0.9863525
## [343,] 0.9758919
## [344,] 0.9807208
## [345,] 0.9698064
## [346,] 0.9908304
## [347,] 0.9935368
## [348,] 0.9767039
## [349,] 0.9892676
## [350,] 0.9768757
## [351,] 0.9697371
## [352,] 0.9830097
## [353,] 0.9816032
## [354,] 0.9796852
## [355,] 0.9671425
## [356,] 0.9609578
## [357,] 0.9882328
## [358,] 0.9909082
## [359,] 0.9691078
## [360,] 0.9873181
## [361,] 0.9759762
## [362,] 0.9910247
## [363,] 0.9975017
## [364,] 0.9835700
## [365,] 0.9900852
## [366,] 0.9876827
## [367,] 0.9773260
## [368,] 0.9866744
## [369,] 0.9880070
## [370,] 0.9868954
## [371,] 0.9875619
## [372,] 0.9852823
## [373,] 0.9787238
## [374,] 0.9834389
## [375,] 0.9849624
## [376,] 0.9891104
## [377,] 0.9899460
## [378,] 0.9805278
## [379,] 0.9867650
## [380,] 0.9865149
## [381,] 0.9737775
## [382,] 0.9877083
## [383,] 0.9779478
## [384,] 0.9915640
## [385,] 0.9815643
## [386,] 0.9783814
## [387,] 0.9877523
## [388,] 0.9831595
## [389,] 0.9766392
## [390,] 0.9910211
## [391,] 0.9821164
## [392,] 0.9765502
## [393,] 0.9828912
## [394,] 0.9816062
## [395,] 0.9870665
## [396,] 0.9826711
## [397,] 0.9885199
## [398,] 0.9853577
## [399,] 0.9848207
## [400,] 0.9810333
## [401,] 0.9905143
## [402,] 0.9857274
## [403,] 0.9791583
## [404,] 0.9885718
## [405,] 0.9824688
## [406,] 0.9850151
## [407,] 0.9697714
## [408,] 0.9895567
## [409,] 0.9799941
## [410,] 0.9829924
## [411,] 0.9799243
## [412,] 0.9720235
## [413,] 0.9805245
## [414,] 0.9925075
## [415,] 0.9763528
## [416,] 0.9755665
## [417,] 0.9810014
## [418,] 0.9942708
## [419,] 0.9800141
## [420,] 0.9854262
## [421,] 0.9780999
## [422,] 0.9924488
## [423,] 0.9848296
## [424,] 0.9711864
## [425,] 0.9803449
## [426,] 0.9862236
## [427,] 0.9833709
## [428,] 0.9877887
## [429,] 0.9914505
## [430,] 0.9810950
## [431,] 0.9925582
## [432,] 0.9881084
## [433,] 0.9875776
## [434,] 0.9943627
## [435,] 0.9826230
## [436,] 0.9761026
## [437,] 0.9860950
## [438,] 0.9814280
## [439,] 0.9858768
## [440,] 0.9912661
## [441,] 0.9914903
## [442,] 0.9885328
## [443,] 0.9867970
## [444,] 0.9823511
## [445,] 0.9897678
## [446,] 0.9713171
## [447,] 0.9898756
## [448,] 0.9661684
## [449,] 0.9732891
## [450,] 0.9864919
## [451,] 0.9707419
## [452,] 0.9877775
## [453,] 0.9620953
## [454,] 0.9892364
## [455,] 0.9892378
## [456,] 0.9865116
## [457,] 0.9921685
## [458,] 0.9715494
## [459,] 0.9931353
## [460,] 0.9834879
## [461,] 0.9790500
## [462,] 0.9909603
## [463,] 0.9889595
## [464,] 0.9866710
## [465,] 0.9875230
## [466,] 0.9760536
## [467,] 0.9794079
## [468,] 0.9935134
## [469,] 0.9797826
## [470,] 0.9831767
## [471,] 0.9767267
## [472,] 0.9900313
## [473,] 0.9963928
## [474,] 0.9854692
## [475,] 0.9821503
## [476,] 0.9751868
## [477,] 0.9866722
## [478,] 0.9867970
## [479,] 0.9901449
## [480,] 0.9902001
## [481,] 0.9706352
## [482,] 0.9813136
## [483,] 0.9861579
## [484,] 0.9826809
## [485,] 0.9897867
## [486,] 0.9887195
## [487,] 0.9647443
## [488,] 0.9653948
## [489,] 0.9843863
## [490,] 0.9806394
## [491,] 0.9824895
## [492,] 0.9816991
## [493,] 0.9818252
## [494,] 0.9825702
## [495,] 0.9815916
## [496,] 0.9885905
## [497,] 0.9922883
## [498,] 0.9808455
## [499,] 0.9756863
## [500,] 0.9833595hist(b3$t)