# NOT RUN {
# The following example follows the example provided at:
#
# http://influentialpoints.com/Training/bootstrap_confidence_intervals.htm
#
# which is provided with a creative commons license:
#
# https://creativecommons.org/licenses/by/3.0/
#
y <- c( 7, 7, 6, 9, 8, 7, 8, 7, 7, 7, 6, 6, 6, 8, 7, 7, 7, 7, 6, 7,
8, 7, 7, 6, 8, 7, 8, 7, 8, 7, 7, 7, 5, 7, 7, 7, 6, 7, 8, 7, 7,
8, 6, 9, 7, 14, 12, 10, 13, 15 )
trm <- function( data, ... ) {
res <- try( mean( data, trim = 0.1, ... ) )
if( class( res ) == "try-error" ) return( NA )
else return( res )
} # end of 'trm' function.
genf <- function( data, par, n, ... ) {
y <- data * par
h <- 1.06 * sd( y ) / ( n^( 1 / 5 ) )
y <- y + rnorm( rnorm( n, 0, h ) )
y <- round( y * ( y > 0 ) )
return( y )
} # end of 'genf' function.
look <- tibber( x = y, statistic = trm, B = 500, rmodel = genf,
test.pars = seq( 0.85, 1.15, length.out = 100 ) )
look
plot( look )
# outer vertical blue lines should cross horizontal blue lines
# near where an estimated p-value is located.
tibber( x = y, statistic = trm, B = 500, rmodel = genf, test.pars = 1 )
# }
# NOT RUN {
look2 <- tibberRM(x = y, statistic = trm, B = 500, rmodel = genf, startval = 1,
step.size = 0.03, verbose = TRUE )
look2
# lower achieved est. p-value should be close to 0.025
# upper should be close to 0.975.
plot( look2 )
trm2 <- function( data, par, n, ... ) {
a <- list( ... )
res <- try( mean( data, trim = a$trim ) )
if( class( res ) == "try-error" ) return( NA )
else return( res )
} # end of 'trm2' function.
tibber( x = y, statistic = trm2, B = 500, rmodel = genf,
test.pars = seq( 0.85, 1.15, length.out = 100 ), trim = 0.1 )
# Try getting the STIB interval. v.terms = 2 below because mfun
# returns the variance of the estimated parameter in the 2nd position.
#
# Note: the STIB interval can be a bit unstable.
mfun <- function( data, ... ) return( c( mean( data ), var( data ) ) )
gennorm <- function( data, par, n, ... ) {
return( rnorm( n = n, mean = mean( data ), sd = sqrt( par ) ) )
} # end of 'gennorm' function.
set.seed( 1544 )
z <- rnorm( 50 )
mean( z )
var( z )
# Trial-and-error is necessary to get a good result with interpolation method.
res <- tibber( x = z, statistic = mfun, B = 500, rmodel = gennorm,
test.pars = seq( 0.95, 1.10, length.out = 100 ), v.terms = 2 )
res
plot( res )
# Much trial-and-error is necessary to get a good result with RM method.
# If it fails to converge, try increasing the tolerance.
res2 <- tibberRM( x = z, statistic = mfun, B = 500, rmodel = gennorm,
startval = c( 0.95, 1.1 ), step.size = 0.003, tol = 0.001, v.terms = 2,
verbose = TRUE )
# Note that it only gives the STIB interval.
res2
plot( res2 )
# }
Run the code above in your browser using DataLab