Hi!
Here is a better code, using the percentil value instead the position, and
some corrections on the graph code.
The reason in the difference on the tails is elementar, but I don't had see
it before. The right and left limits are calculated by simulation, and his
differences from the mean are not exactly equal, so the areas under the
normal curve! I'm feeling like dummie now...
Thanks for the help!
Raphael Saldanha
BRAZIL
On Wed, Aug 20, 2008 at 12:07 PM, Raphael Saldanha <
[EMAIL PROTECTED]> wrote:
> 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."
>
--
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.
