Thanks Professor and Bernardo,
What I'm trying to do is this: I have a macro for Minitab. His author says
it's a "Monte Carlo simulation to estimate a confidence interval". But and
don't have Minitab, don't like to work with illegal licenses, and LOVE R.
So, I re-write the macro in a R script.
Bellow, is the original macro, for Minitab:
GMACRO
TesteMedia
do K1=1:1000
random 200 C1;
normal 100 2.
mean C1 K2
let C2(K1)=K2
enddo
sort C2 C3
let K3=C3(25)
let K4=C3(975)
print K3 K4
ENDMACRO
The macro's author think on this way: I have a y vector of sorted means with
1000 registers. So, my CI is between the 25º and 975º element.
So, I ask: Is this a Monte Carlo Simulation, and the nome is correct? The
way to isolate the inferior and superior means is correct?
About the graphics, I know the sample can't be reproduced because it's
random. But, re-running the script some times, I can see, some times,
differences in the right and left tail, under the normal curve. So, I wonder
to know where my code is wrong, but just for know, I agree with histogram
isn't the best way to illustrate confidence interval. I just show him on the
plot to illustrated the sample used. I will re-make the plot, just with the
curve and areas after the confidence interval, without the histogram.
Thanks for the help!
On Wed, Aug 20, 2008 at 12:16 AM, Prof Brian Ripley
<[EMAIL PROTECTED]>wrote:
> 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/<http://www.stats.ox.ac.uk/%7Eripley/>
> University of Oxford, Tel: +44 1865 272861 (self)
> 1 South Parks Road, +44 1865 272866 (PA)
> Oxford OX1 3TG, UK Fax: +44 1865 272595
>
--
Atenciosamente,
Raphael Saldanha
[EMAIL PROTECTED]
Robert Orben - "To err is human - and to blame it on a computer is even
more so."
[[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.
