On Tue, 19 Aug 2008, Raphael Saldanha wrote:
Hi!
With the following script, I'm trying to make a demonstration of a
Confidence Interval, but I'm observing some differences on tails.
You need to tell us what you are trying to do. You seem to be computing a
parametric bootstrap interval, but incorrectly (you need to reflect
percentile intervals). See Davison & Hinkley (1997) 'The Bootstrap and
its Application' for more details.
In any case, your simulation is not repeatable (you set no seed), so we
don't know what you saw and what 'differences' disturbed you.
Your calculation of quantiles is not quite correct: use quantile().
(The indices go from 1 to 1000, not 0 to 1000.) E.g.
quantile(1:1000, c(0.025, 0.975))
2.5% 97.5%
25.975 975.025
and not 25, 975.
As your results show, a histogram is not a good summary of the shape of a
distribution on 1000 points. We can do much better with an ecdf or a
density estimate.
# Teste de m?dia entre uma amostra e uma popula??o normal
# Autor: Raphael de Freitas Saldanha
# Agosto de 2008
n <- 200 # Sample size
xbar <- 100 # Sample mean
s <- 2 # Sample SD
nc <- 0.95 # Confidence level (95% -> 0.95)
rep <- 1000 # Loops
#######################################################################
y <- NULL # Vetor com as m?dias da amostra
for (i in 1:rep){ # Loop
x <- rnorm(n,xbar,s) # Gere uma amostra normal n elementos, xbar m?dia e
s desvio-padr?o
x <- mean(x) # Calcule a m?dia (exata) dessa amostra
y <- c(y,x) # Coloque essa m?dia em um registro em y
}
y <- sort(y) # Ordene todas as m?dias geradas
LI <- y[((1-nc)/2)*rep] # Limite inferior,
LS <- y[rep-(((1-nc)/2)*rep)] # Limite superior
### PARTE GR?FICA ###
x <- mean(y)
xvals <- seq(-LI, LI, length.out=5000)
dvals <- dnorm(xvals,mean(y), sd(y))[1:5000]
xbvals <- seq(LS, LS*2, length.out=5000)
dbvals <- dnorm(xbvals,mean(y), sd(y))[1:5000]
ahist <- hist(y, freq=FALSE, col="lightblue", main="Intervalo de confian?a")
polygon(c(xvals,LI,LI), c(dvals,dnorm(LI,mean(y), sd(y)),min(dvals)),
col="orange", lty=0)
polygon(c(LS,LS,xbvals), c(min(dbvals),dnorm(LS,mean(y), sd(y)),dbvals),
col="orange", lty=0)
curve(dnorm(x,mean(y), sd(y)),add=TRUE, lty=1, col="red",lwd=2)
### Intervalo de Confian?a ###
LI # Limite inferior
LS # Limite superior
--
Atenciosamente,
Raphael Saldanha
[EMAIL PROTECTED]
[[alternative HTML version deleted]]
--
Brian D. Ripley, [EMAIL PROTECTED]
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@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.
