[R] adding one line to a plot

[Dimitri Liakhovitski]

Thanks a lot, Peter, that's exactly what I was looking for:

plot(density(tstat),xlim=c(-5,8),ylim=c(0,.4),lwd=2,col='red')
z <- tstat[order(tstat)]
lines(z,dt(z,df=18),col='blue')
legend(4,.3,c("exact","t(18)"),lwd=c(2,1),col=c('red','blue'))

Dimitri


On Mon, May 24, 2010 at 1:26 PM, Peter Ehlers <ehlers at ucalgary.ca> wrote:
> On 2010-05-24 10:50, Dimitri Liakhovitski wrote:
>>
>> Hello!
>>
>> I am running a very simple mini Monte-Carlo below using the function
>> tstatistic (right below this sentence):
>>
>> tstatistic = function(x,y){
>>        m=length(x)
>>        n=length(y)
>>        sp=sqrt( ((m-1)*sd(x)^2 + (n-1)*sd(y)^2)/(m+n-2) )
>>        t=(mean(x)-mean(y))/(sp*sqrt(1/m+1/n))
>>        return(t)
>> }
>>
>> alpha=.1; m=10; n=10 # sets alpha, m, n - for run 1
>> N=10000 # sets the number of simulations
>> n.reject=0 # counter of num. of rejections
>> tstat<-vector()
>> for (i in 1:N){
>>        x = rnorm(m,mean=10,sd=2) # simulates xs from population 1
>>        y=rexp(n,rate=1/10)
>>        t = tstatistic(x,y) # computes the t statistic
>>        tstat = c(tstat,t)
>>        if (abs(t)>qt(1-alpha/2,n+m-2))
>>        n.reject=n.reject+1 # reject if |t| exceeds critical pt
>> }
>> true.sig.level=n.reject/N # est. is proportion of rejections
>> true.sig.level
>>
>> And then I would like to plot the observed t values (in the vector
>> "tstat") against the values of the t density with df=18. Somehow, my t
>> denstity (code line 2 below) does not show up:
>> plot(density(tstat),xlim=c(-5,8),ylim=c(0,.4),lwd=3)
>> lines(x,dt(x,df=18))
>> legend(4,.3,c("exact","t(18)"),lwd=c(3,1))
>>
>
> Perhaps
>
>  plot(x,dt(x,df=18))
>
> will provide a clue. (What's x?)
> You probably want
>
>  x <- tstat[order(tstat)]
>
>  -Peter Ehlers
>



-- 
Dimitri Liakhovitski
Ninah Consulting
www.ninah.com