[R] adding one line to a plot

Peter Ehlers ehlers at ucalgary.ca
Mon May 24 19:26:09 CEST 2010


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



More information about the R-help mailing list