In order to apply the bootstrap, you must resample, uniformly at random
from the independent units of measurement in your data. Assuming that
these represent the rows of 'data', consider the following:
est <- function(y, x, obeta = c(1,1), verbose=FALSE) {
n <- length(x)
X <- cbind(rep(1, n), x)
nbeta <- c(0,0)
iter <- 0
while(crossprod(obeta-nbeta)>10^(-12)) {
nbeta <- obeta
eta <- X%*%nbeta
mu <- eta
mu1 <- 1/eta
W <- diag(as.vector(mu1))
Z <- X%*%nbeta+(y-mu)
XWX <- t(X)%*%W%*%X
XWZ <- t(X)%*%W%*%Z
Cov <- solve(XWX)
obeta <- Cov%*%XWZ
iter <- iter+1
if(verbose)
cat("Iteration # and beta1= ",iter, nbeta, "\n")
}
return(nbeta[1,1])
}
boot <- function(data, reps) {
n <- nrow(data)
Nt <- vector('numeric', length=reps)
for(Ncount in 1:reps) {
#resample the rows of data
bdata <- data[sample(1:n,n,replace=TRUE),]
#recompute and store estimate
Nt[Ncount] <- est(bdata[,1], bdata[,2])
}
return(Nt)
}
stem(boot(data,1000),width=60)
The decimal point is at the |
-3 | 4
-2 |
-1 | 2
-0 | 888666666555444333333333332222222222211111111111
0 | 000000000000001111111111111111111111111111222222+400
1 | 000000000000000000000000000000011111111111111111+203
2 | 000000000000000011111111222222222223333444444444+23
3 | 000000001111111111222222222223344444445555556666
4 | 11222233445556666789
5 | 02334446677899
6 | 1112334455778
7 | 000011235568
8 | 001799
9 | 0259
10 | 1446
11 | 19
12 | 48
13 | 8
14 | 024
15 |
16 |
17 | 0788
18 |
19 | 1
On Wed, 2011-07-20 at 18:09 -0400, Val wrote:
> Hi all,
>
> I am facing difficulty on how to use bootstrap sampling and
> below is my example of function.
>
> Read a data , use some functions and use iteration to find the solution(
> ie, convergence is reached). I want to use bootstrap approach to do it
> several times (200 or 300 times) this whole process and see the
> distribution of parameter of interest.
>
> Below is a small example that resembles my problem. However, I found out
> all samples are the same. So I would appreciate your help on this case.
>
> #**************************************
> rm(list=ls())
> xx <- read.table(textConnection(" y x
> 11 5.16
> 11 4.04
> 14 3.85
> 19 5.68
> 4 1.26
> 23 7.89
> 15 4.25
> 17 3.94
> 7 2.35
> 17 4.74
> 14 5.49
> 11 4.12
> 17 5.92"), header=TRUE)
> data <- as.matrix(xx)
> closeAllconnections()
>
> Nt <- NULL
> for (Ncount in 1:100)
> {
> y <- data[,1]
> x <- data[,2]
> n <- length(x)
>
> X <- cbind(rep(1,n),x) #covariate/design matrix
> obeta<- c(1,1) #previous/starting values of beta
>
> nbeta <- c(0,0) #new beta
> iter=0
>
> while(crossprod(obeta-nbeta)>10^(-12))
> {
> nbeta <- obeta
> eta <- X%*%nbeta
> mu <- eta
> mu1 <- 1/eta
> W <- diag(as.vector(mu1))
> Z <- X%*%nbeta+(y-mu)
> XWX <- t(X)%*%W%*%X
> XWZ <- t(X)%*%W%*%Z
> Cov <- solve(XWX)
> obeta <- Cov%*%XWZ
> iter <- iter+1
>
> cat("Iteration # and beta1= ",iter, nbeta, "\n")
> }
>
> Nt[Ncount] <- nbeta[1,1]
> }
> Nt
> summary(Nt)
> #**************e*****************************************
>
> [[alternative HTML version deleted]]
>
> ______________________________________________
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
