[R] Bootstrap Confidence Intervals
Chuck Cleland
ccleland at optonline.net
Tue Jan 1 15:39:01 CET 2008
_Fede_ wrote:
> Hi again.
>
> Watching this example that appears in the help page
>
> ratio <- function(d, w) sum(d$x * w)/sum(d$u * w)
> city.boot <- boot(city, ratio, R = 999, stype = "w",sim = "ordinary")
> boot.ci(city.boot, conf = c(0.90,0.95),type =
> c("norm","basic","perc","bca"))
>
> I have tried to do the following (calling boot() to create an object to pass
> to boot.ci):
>
> x <- rnorm(20)
> kurtosis <- function(x) (mean((x-mean(x))^4))/(sd(x)^4)
> x.boot <- boot(x, kurtosis, R = 999, sim = "ordinary")
> boot.ci(x.boot, conf = 0.95,type = c("norm","basic","perc","bca"))
>
> But I don't know why this don't work. The editor window shows the following
> error message:
>
> Error in statistic(data, original, ...) : unused argument(s) (1:20)
>
> I suppose that something is wrong with my data but I don't know what is.
>
> Thanks in advance and wishing everybody a happy new year.
Check the statistic argument to boot() very carefully. When sim =
"ordinary" the statistic function must have at least two arguments. Try
something like this:
library(boot)
kurtosis <- function(x) (mean((x-mean(x))^4))/(sd(x)^4)
x <- rnorm(20)
x.boot <- boot(x,
statistic = function(d, ind){kurtosis(d[ind])},
R = 999,
sim = "ordinary")
boot.ci(x.boot, conf = 0.95, type = c("norm","basic","perc","bca"))
BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 999 bootstrap replicates
CALL :
boot.ci(boot.out = x.boot, conf = 0.95, type = c("norm", "basic",
"perc", "bca"))
Intervals :
Level Normal Basic
95% ( 1.060, 2.430 ) ( 0.899, 2.233 )
Level Percentile BCa
95% ( 1.394, 2.728 ) ( 1.373, 2.690 )
Calculations and Intervals on Original Scale
Also, note that the e1071 package contains a kurtosis function.
> Regards
>
> _Fede_
>
> Prof Brian Ripley wrote:
>> You need to call boot() to create an object to pass to boot.ci().
>>
>> There are lots of examples in the help pages and in the book that package
>> 'boot' supports. From the help:
>>
>> Usage:
>>
>> boot.ci(boot.out, conf = 0.95, type = "all",
>> index = 1:min(2,length(boot.out$t0)), var.t0 = NULL,
>> var.t = NULL, t0 = NULL, t = NULL, L = NULL, h = function(t)
>> t,
>> hdot = function(t) rep(1,length(t)), hinv = function(t) t,
>> ...)
>>
>> Arguments:
>>
>> boot.out: An object of class '"boot"' containing the output of a
>> bootstrap calculation.
>>
>> and try class(z) .
>>
>>
>> On Sun, 30 Dec 2007, _Fede_ wrote:
>>
>>>
>>> Hi all.
>>>
>>> This is my first post in this forum. Finally I find a forum in the web
>>> about
>>> R, although is not in my language.
>>>
>>> Now I'm working with Bootstrap CI. I'd like to know how I can calculate a
>>> Bootstrap CI for any statistic, in particular, for Kurtosis Coeficient. I
>>> have done the following code lines:
>>>
>>>> library(boot)
>>>> x=rnorm(20)
>>>> kurtosis=function(x) (mean((x-mean(x))^4))/(sd(x)^4)
>>>> z <- numeric(10000)
>>>> for(i in 1:10000)
>>>> z[i]=kurtosis(sample(x, replace=TRUE))
>>>> boot.ci(z, conf = 0.95,type = c("norm","basic","perc","bca"))
>>> But the output shows the next error:
>>>
>>> Error en if (ncol(boot.out$t) < max(index)) { :
>>> argumento tiene longitud cero
>>>
>>> I don't know what is wrong.
>>>
>>> I hope that somebody can help me. Sorry for my english.
>>>
>>> All have a nice new year.
>>>
>>> _Fede_
>>>
>> --
>> Brian D. Ripley, ripley at stats.ox.ac.uk
>> Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/
>> University of Oxford, Tel: +44 1865 272861 (self)
>> 1 South Parks Road, +44 1865 272866 (PA)
>> Oxford OX1 3TG, UK Fax: +44 1865 272595
>>
>> ______________________________________________
>> R-help at 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.
--
Chuck Cleland, Ph.D.
NDRI, Inc.
71 West 23rd Street, 8th floor
New York, NY 10010
tel: (212) 845-4495 (Tu, Th)
tel: (732) 512-0171 (M, W, F)
fax: (917) 438-0894
More information about the R-help
mailing list