Hi I'm trying to plot the power functions of a t-test and a sign test using simulated data from a normal distribution N(theta,1). The power function of the t-test is Pr(TS1>c1) and the power function of the sign test is Pr(TS2>c2). where TS1 is the test statistic of the t-test which is mean(x)/(sd(x)*sqrt(n)) and TS2 is the test statistic of the sign test which is sum(x>0). I have to do this for 0<theta<2
I have produced the following code but can't figure out where I'm going wrong. Any help would be appreciated - thanks!
c1 <- 1.9
c2 <- 8
n <- 10
thetaVals = seq(0,2,length.out = 1000)
set.seed(1)
P1 <- vector("numeric", length = 1000)
P2 <- vector("numeric", length = 1000)
estimators <- vector("numeric", length = 1000)
for (n in seq(1,10)) {
for (i in 1:1000) {
x <- rnorm(thetaVals, sd = 1)
TS1[i] = mean(x)/(sd(x)*sqrt(n))
TS2[i] = sum(x>0)
}
P1[n] <- Pr(TS1 > c1)
P2[n] <- PR(TS2 > c2)
}
nVal <- seq(1, 1000 , by = 1)
First thing I noticed is that TS1 and TS2 are not defined. Define them as you did for P1 and P2 before the loop.
I am not familiar with the Pr function, but you have written Pr and PR, this might also cause a problem. If I understand this correctly, you could substitute Pr by mean, what will give you the proportion of TS > c.
Hi, stupid mistake from me.
Pr(TS1 > c1) isn't executable code. It was just me trying to say probability that TS1 exceeds c1.
Any way I can convert this into code?
So this is my updated code but for some reason when I plot Power (P1) against the theta values i'm not getting the correct plot, instead i'm getting a flat line.
c1 <- 1.5
c2 <- 6
n <- 10
set.seed(1) #Ensures same random numbers are generated each time
P1 <- vector("numeric", length = 1000) #Create vector for first power function
P2 <- vector("numeric", length = 1000) #Create vector for second power function
TS1 <- vector("numeric", length = 1000) #Create vector for first test statistic
TS2 <- vector("numeric", length = 1000) #Create vector for second test statistic
for (theta in seq(0,2)) {
for (i in 1:1000) {
y <- rnorm(1000, theta, 1) #Generate random values of y from normal distribution using mean = theta, sd = 1
TS1[i] = mean(y)*sqrt(n)/(sd(y)) #Define first test statistic
TS2[i] = sum(y > 0) #Define second test statistic
}
P1[theta] <- mean(TS1 > c1) #Calculates power for each test statistic for power function 1
P2[theta] <- mean(TS2 > c2) #Calculates power for each test statistic for power function 2
}
thetaVals <- seq(0, 2, length.out =1000)
plot(thetaVals, P1, type ="p", xlab ="theta")
plot(thetaVals, P2, type ="p", xlab ="theta")