Remuestreo - Bootstrap
Eduardo Lezama Meneses
3/9/2020
R Markdown
This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.
When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:
summary(cars)## speed dist
## Min. : 4.0 Min. : 2.00
## 1st Qu.:12.0 1st Qu.: 26.00
## Median :15.0 Median : 36.00
## Mean :15.4 Mean : 42.98
## 3rd Qu.:19.0 3rd Qu.: 56.00
## Max. :25.0 Max. :120.00Including Plots
You can also embed plots, for example:
set.seed(124)
MO<- sort(rnorm(120,2,0.2))
CIC<- sort(rnorm(120,10,1.5))
plot(MO,CIC) 
# Grafico descriptivocor(MO, CIC)## [1] 0.9944214## Ajustando un modelo de regresión simple
model<-lm(CIC~MO);summary(model)##
## Call:
## lm(formula = CIC ~ MO)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.53490 -0.06721 -0.00820 0.06231 0.89548
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -8.05203 0.17600 -45.75 <2e-16 ***
## MO 8.98922 0.08778 102.41 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1635 on 118 degrees of freedom
## Multiple R-squared: 0.9889, Adjusted R-squared: 0.9888
## F-statistic: 1.049e+04 on 1 and 118 DF, p-value: < 2.2e-16\[\hat {CIC} = -8.05 + 8.99* MO\]
## Prueba de correlación
\[H_0: \rho = 0 \\ H_a: \rho \neq 0\]
pr1<-cor.test(CIC, MO, alternative = "t",method = "pearson"); pr1##
## Pearson's product-moment correlation
##
## data: CIC and MO
## t = 102.41, df = 118, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.9919946 0.9961139
## sample estimates:
## cor
## 0.9944214ifelse ( pr1$p.value <0.05, "Rechazo H0", "No rechazo H0")## [1] "Rechazo H0"### Entendiendo Bootstrap
library(gtools)
data=c(260,195,335)
#todas las posibles muestras de tamaño 3 tal como el vector
df = data.frame(permutations(n = 3, r = 3, v = data, repeats.allowed = T))
View(df)
rowMeans(df)## [1] 195.0000 216.6667 241.6667 216.6667 238.3333 263.3333 241.6667 263.3333
## [9] 288.3333 216.6667 238.3333 263.3333 238.3333 260.0000 285.0000 263.3333
## [17] 285.0000 310.0000 241.6667 263.3333 288.3333 263.3333 285.0000 310.0000
## [25] 288.3333 310.0000 335.0000df2 = data.frame(df, rowMeans(df))
colnames(df2) = c('x1*','x2*','x3*','media(x*)')
View(df2)
hist(df2$'media(x*)', xlim = c(100,400))# distribución bootstrap
rug(df2$'media(x*)')
abline(v = mean(df2$'media(x*)'), col = 'red', lwd = 2, lty = 2)
## Error estandar de las medias
eeb =sd(df2$'media(x*)')
eeb## [1] 33.65549media = mean(df2$'media(x*)')
hist(df2$'media(x*)', xlim = c(100,400))# distribución bootstrap
rug(df2$'media(x*)')
abline(v = mean(df2$'media(x*)'), col = 'red', lwd = 2, lty = 2)
abline(v = media + eeb * c(-1, 1), col = 'blue', lwd = 2, lty = 2)
abline(v = media + eeb * c(-2, 2), col = 'darkgreen', lwd = 2, lty = 2)
eeb =sd(df2$'media(x*)')
eeb## [1] 33.65549hist(df2$'media(x*)', xlim = c(100,400))# distribución bootstrap
rug(df2$'media(x*)', col = "darkgreen", lwd = 2)
abline(v = mean(df2$'media(x*)'), col = 'red', lwd = 2, lty = 2)
abline(v = mean(df2$'media(X*)')+eeb*c(-1,1),col="blue",lwd=2,lty=2)## Warning in mean.default(df2$"media(X*)"): argument is not numeric or logical:
## returning NAabline(v = mean(df2$'media(X*)')+eeb*c(-2,2),col="gold",lwd=2,lty=2)## Warning in mean.default(df2$"media(X*)"): argument is not numeric or logical:
## returning NAsegments(x0 = mean(df2$'media (X*)')-2*eeb,y0 = 5,x1 =mean(df2$'Medias (X*)')+2*eeb, y1 = 5 )## Warning in mean.default(df2$"media (X*)"): argument is not numeric or logical:
## returning NA## Warning in mean.default(df2$"Medias (X*)"): argument is not numeric or logical:
## returning NAtext(x = 210, y= 5.5, "2 EE")
Distribución bootstrap
library(gtools)
data=c(260,195,335, 320)
#todas las posibles muestras de tamaño 4 tal como el vector
df = data.frame(permutations(n = 4, r = 4, v = data, repeats.allowed = T))
View(df)
rowMeans(df)## [1] 195.00 211.25 226.25 230.00 211.25 227.50 242.50 246.25 226.25 242.50
## [11] 257.50 261.25 230.00 246.25 261.25 265.00 211.25 227.50 242.50 246.25
## [21] 227.50 243.75 258.75 262.50 242.50 258.75 273.75 277.50 246.25 262.50
## [31] 277.50 281.25 226.25 242.50 257.50 261.25 242.50 258.75 273.75 277.50
## [41] 257.50 273.75 288.75 292.50 261.25 277.50 292.50 296.25 230.00 246.25
## [51] 261.25 265.00 246.25 262.50 277.50 281.25 261.25 277.50 292.50 296.25
## [61] 265.00 281.25 296.25 300.00 211.25 227.50 242.50 246.25 227.50 243.75
## [71] 258.75 262.50 242.50 258.75 273.75 277.50 246.25 262.50 277.50 281.25
## [81] 227.50 243.75 258.75 262.50 243.75 260.00 275.00 278.75 258.75 275.00
## [91] 290.00 293.75 262.50 278.75 293.75 297.50 242.50 258.75 273.75 277.50
## [101] 258.75 275.00 290.00 293.75 273.75 290.00 305.00 308.75 277.50 293.75
## [111] 308.75 312.50 246.25 262.50 277.50 281.25 262.50 278.75 293.75 297.50
## [121] 277.50 293.75 308.75 312.50 281.25 297.50 312.50 316.25 226.25 242.50
## [131] 257.50 261.25 242.50 258.75 273.75 277.50 257.50 273.75 288.75 292.50
## [141] 261.25 277.50 292.50 296.25 242.50 258.75 273.75 277.50 258.75 275.00
## [151] 290.00 293.75 273.75 290.00 305.00 308.75 277.50 293.75 308.75 312.50
## [161] 257.50 273.75 288.75 292.50 273.75 290.00 305.00 308.75 288.75 305.00
## [171] 320.00 323.75 292.50 308.75 323.75 327.50 261.25 277.50 292.50 296.25
## [181] 277.50 293.75 308.75 312.50 292.50 308.75 323.75 327.50 296.25 312.50
## [191] 327.50 331.25 230.00 246.25 261.25 265.00 246.25 262.50 277.50 281.25
## [201] 261.25 277.50 292.50 296.25 265.00 281.25 296.25 300.00 246.25 262.50
## [211] 277.50 281.25 262.50 278.75 293.75 297.50 277.50 293.75 308.75 312.50
## [221] 281.25 297.50 312.50 316.25 261.25 277.50 292.50 296.25 277.50 293.75
## [231] 308.75 312.50 292.50 308.75 323.75 327.50 296.25 312.50 327.50 331.25
## [241] 265.00 281.25 296.25 300.00 281.25 297.50 312.50 316.25 296.25 312.50
## [251] 327.50 331.25 300.00 316.25 331.25 335.00df2 = data.frame(df, rowMeans(df))
colnames(df2) = c('x1*','x2*','x3*','x4*','media(x*)')
View(df2)
hist(df2$'media(x*)', xlim = c(100,400))# distribución bootstrap
rug(df2$'media(x*)')
abline(v = mean(df2$'media(x*)'), col = 'red', lwd = 2, lty = 2)
eeb =sd(df2$'media(x*)')
eeb## [1] 27.69583media = mean(df2$'media(x*)')
hist(df2$'media(x*)', xlim = c(100,400))# distribución bootstrap
rug(df2$'media(x*)', col = "darkgreen", lwd = 2)
abline(v = mean(df2$'media(x*)'), col = 'red', lwd = 2, lty = 2)
abline(v = mean(df2$'media(X*)')+eeb*c(-1,1),col="blue",lwd=2,lty=2)## Warning in mean.default(df2$"media(X*)"): argument is not numeric or logical:
## returning NAabline(v = mean(df2$'media(X*)')+eeb*c(-2,2),col="gold",lwd=2,lty=2)## Warning in mean.default(df2$"media(X*)"): argument is not numeric or logical:
## returning NAsegments(x0 = mean(df2$'media (X*)')-2*eeb,y0 = 5,x1 =mean(df2$'Medias (X*)')+2*eeb, y1 = 5 )## Warning in mean.default(df2$"media (X*)"): argument is not numeric or logical:
## returning NA## Warning in mean.default(df2$"Medias (X*)"): argument is not numeric or logical:
## returning NAtext(x = 210, y= 5.5, "2 EE")
set.seed(1964)
MO = sort(rnorm(90, 2, 0.25))
pH = sort(runif(90, 5.2, 6.2))
quimi = data.frame(MO, pH)
cor(MO, pH)## [1] 0.9460442fcorr <- function(datos, obs){
d2 <- quimi[obs, ]
return(cor(d2$MO, d2$pH))
}library(boot)##
## Attaching package: 'boot'## The following objects are masked from 'package:gtools':
##
## inv.logit, logitset.seed(1964)
bootcorrP <- boot(quimi, fcorr, R = 500)
bootcorrP##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = quimi, statistic = fcorr, R = 500)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 0.9460442 0.001354613 0.007254421mean(bootcorrP$t)## [1] 0.9473989sd(bootcorrP$t)## [1] 0.007254421plot(bootcorrP)
boot.ci(bootcorrP, type = c('norm', 'basic', 'perc', 'bca'))## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 500 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = bootcorrP, type = c("norm", "basic", "perc",
## "bca"))
##
## Intervals :
## Level Normal Basic
## 95% ( 0.9305, 0.9589 ) ( 0.9312, 0.9592 )
##
## Level Percentile BCa
## 95% ( 0.9329, 0.9609 ) ( 0.9224, 0.9571 )
## Calculations and Intervals on Original Scale
## Some BCa intervals may be unstableset.seed(1964)
C = sort(rnorm(90, 0.15, 0.02))
N = sort(rnorm(90, 0.30, 0.05))
rCN=C/N
quimi = data.frame(C, N, rCN)
head(quimi)(media = mean(quimi$rCN))## [1] 0.4948643fcorr <- function(datos, obs){
d2 <- quimi$rCN[obs, ]
return(cor(d2$MO, d2$pH))
}library(boot)
set.seed(1964)
MO = sort(rnorm(90, 2, 0.25))
pH = sort(runif(90, 5.2, 6.2))
quimi = data.frame(MO, pH)
cor(MO, pH)## [1] 0.9460442fcorr <- function(datos, obs){
d2 <- quimi[obs, ]
return(cor(d2$MO, d2$pH))
}
library(boot)
#
set.seed(1964)
bootcorrP <- boot(quimi, fcorr, R = 500)
bootcorrP##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = quimi, statistic = fcorr, R = 500)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 0.9460442 0.001354613 0.007254421mean(bootcorrP$t)## [1] 0.9473989sd(bootcorrP$t)## [1] 0.007254421plot(bootcorrP)
#library(growthmodels)
growth <- chapmanRichards(0:10, 10, 0.5, 0.3, 0.5)
# Calculate inverse function
time <- chapmanRichards.inverse(growth, 10, 0.5, 0.3, 0.5)
plot(growth, pch=10)
definimos la función
Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.