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))
Thanks a lot!
D.
adding one line to a plot
4 messages · Peter Ehlers, Dimitri Liakhovitski, Jannis
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
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
1 day later
Your line does show up, but the x values start at ~6 and the y values are incredibly small, so it is masked by the y=0 line. You probably have to check your calculations! --- Dimitri Liakhovitski <dimitri.liakhovitski at gmail.com> schrieb am Mo, 24.5.2010:
Von: Dimitri Liakhovitski <dimitri.liakhovitski at gmail.com>
Betreff: [R] adding one line to a plot
An: r-help at r-project.org
Datum: Montag, 24. Mai, 2010 16:50 Uhr
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))
Thanks a lot!
D.
______________________________________________ R-help at 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.