R/functions.R
In StrattonCh/powerapppackage:
exp_samp_max_pwrfunc_greater <- function(theta, alpha, theta_not, n){
#function to calculate power function of sample max from exponential(theta) distribution
1 - (1 - exp(theta_not*log(1 - (1 - alpha)^(1/n)) / (theta) ))^n
}
exp_samp_max_pwrfunc_less <- function(theta, alpha, theta_not, n){
#function to calculate power function of sample max from exponential(theta) distribution
(1 - exp(theta_not*log(1 - (alpha)^(1/n)) / (theta) ))^n
}
exp_samp_max_pwrfunc_noteqto <- function(theta, alpha, theta_not, n){
#function to calculate power function of sample max from exponential(theta) distribution
1 - (1 - exp(theta_not*log(1 - (1 - alpha/2)^(1/n)) / (theta) ))^n + (1 - exp(theta_not*log(1 - (alpha/2)^(1/n)) / (theta) ))^n
}
norm_samp_min_pwrfunc_greater <- function(theta, alpha, sigma, theta_not, n){
#function to calculate power of sample min of normal(theta, sigma^2) distribution
k <- qnorm(1 - alpha^(1/n), theta_not, sigma)
(1 - pnorm(k, theta, sigma))^n
}
norm_samp_min_pwrfunc_less <- function(theta, alpha, sigma, theta_not, n){
#function to calculate power of sample min of normal(theta, sigma^2) distribution
k <- qnorm(1 - (1 - alpha)^(1/n), theta_not, sigma)
1 - (1 - pnorm(k, theta, sigma))^n
}
norm_min_cdf <- function(x, theta, sigma, n){
#function to calculate cdf of sample min from normal(theta, sigma^2)
1 - (1 - pnorm(x, theta, sigma))^n
}
norm_samp_min_pwrfunc_noteqto <- function(theta, alpha, sigma, theta_not, n){
#function to calculate power of sample min of normal(theta, sigma^2) distribution
k1 <- qnorm(1 - (1 - alpha/2)^(1/n), theta_not, sigma)
k2 <- qnorm(1 - (alpha/2)^(1/n), theta_not, sigma)
1 - norm_min_cdf(k2, theta, sigma, n) + norm_min_cdf(k1, theta, sigma, n)
}
norm_max_cdf <- function(x, theta, sigma, n){
#function to calculate the cdf of the sample max from normal(theta, sigma^2)
(pnorm(x, theta, sigma))^n
}
norm_samp_max_pwrfunc_greater <- function(theta, alpha, sigma, theta_not, n){
#function to calculate power of sample max of normal(theta, sigma^2) distribution
k <- qnorm((1 - alpha)^(1/n), theta_not, sigma)
1 - norm_max_cdf(k, theta, sigma, n)
}
norm_samp_max_pwrfunc_less <- function(theta, alpha, sigma, theta_not, n){
#function to calculate power of sample max of normal(theta, sigma^2) distribution
k <- qnorm(alpha^(1/n), theta_not, sigma)
norm_max_cdf(k, theta, sigma, n)
}
norm_samp_max_pwrfunc_noteqto <- function(theta, alpha, sigma, theta_not, n){
#function to calculate power of sample max of normal(theta, sigma^2) distribution
k1 <- qnorm((alpha/2)^(1/n), theta_not, sigma)
k2 <- qnorm((1 -alpha/2)^(1/n), theta_not, sigma)
1 - norm_max_cdf(k2, theta, sigma, n) + norm_max_cdf(k1, theta, sigma, n)
}
unif_samp_max_pwrfunc_greater <- Vectorize(function(theta, alpha, theta_not, n){
#function to calculate power of sample max of unif(0,theta)
k <- theta_not*(1 - alpha)^(1/n)
if(is.na(theta)){
return(NA)
}
if(theta < k){
return(0)
}
if(k <= theta){
return(1 - k^n/theta^n)
}
})
unif_samp_max_pwrfunc_less <- Vectorize(function(theta, alpha, theta_not, n){
#function to calculate power of sample max of unif(0,theta)
k <- theta_not*(alpha)^(1/n)
if(is.na(theta)){
return(NA)
}
if(theta < k){
return(1)
}
if(k <= theta){
return(k^n/theta^n)
}
})
unif_samp_max_pwrfunc_noteqto <- Vectorize(function(theta, alpha, theta_not, n){
#function to calculate power of sample max of unif(0,theta)
k1 <- theta_not*(alpha/2)^(1/n)
k2 <- theta_not*(1 - alpha/2)^(1/n)
if(is.na(theta)){
return(NA)
}
if(theta <= k1){
return(1)
}
if(k1 < theta & theta <= k2){
return(alpha*theta_not^n / (2*theta^n))
}
if(k2 <= theta){
return(1 - theta_not^n * (1 - alpha/2) / theta^n + alpha*theta_not^n / (2*theta^n))
}
})
unif_samp_min_pwrfunc_greater <- Vectorize(function(theta, alpha, theta_not, n){
#function to calculate power of sample min of unif(0,theta)
k <- theta_not*(1 - alpha^(1/n))
if(is.na(theta)){
return(NA)
}
if(theta <= k){
return(0)
}
if(theta > k){
return((1 - k/theta)^n)
}
})
unif_samp_min_pwrfunc_less <- Vectorize(function(theta, alpha, theta_not, n){
#function to calculate power of sample min of unif(0,theta)
k <- theta_not*(1 - (1 - alpha)^(1/n))
if(is.na(theta)){
return(NA)
}
if(theta <= k){
return(1)
}
if(theta > k){
return(1 - (1 - k/theta)^n)
}
})
unif_samp_min_pwrfunc_noteqto <- Vectorize(function(theta, alpha, theta_not, n){
#function to calculate power of sample min of unif(0,theta)
k1 <- theta_not*(1 - (1 - alpha/2)^(1/n))
k2 <- theta_not*(1 - (alpha/2)^(1/n))
if(is.na(theta)){
return(NA)
}
if(theta <= k1){
return(1)
}
if(k1 < theta & theta <= k2){
return(1 - (1 - k1/theta)^n)
}
if(theta > k2){
return(1 + (1 - k2/theta)^n - (1 - k1/theta)^n)
}
})
dirwinhall <- Vectorize(function(x, n, theta){
# function to calculate density of sum of n unif(0, theta) RVs
# inputs : x - value at which to calculate density
# : n - number of unif RVs to sum
# : theta - population maximum
# outputs : numeric value that is the density at x
if(x/theta < 0) return(0)
if(x/theta > n) return(0)
if(x/theta >= 0 & x/theta <= n){
X <- floor(x/theta)
r <- seq(from = 0, to = X)
s <- (-1)^r * choose(n, r)*(x/theta-r)^(n-1)/factorial(n-1)
return(sum(s)/theta)
}
})
pirwinhall <- Vectorize(function(q, n, theta){
# function to calculate cumulative density of sum of n unif(0, theta) RVs
# inputs : q - quantile
# : n - number of unif RVs to sum
# : theta - population maximum
# outputs : numeric value that is the cumulative density at x
if(q/theta < 0) return(0)
if(q/theta > n) return(1)
if(q/theta >= 0 & q/theta <= n){
X <- floor(q/theta)
r <- seq(from = 0, to = X)
s <- (-1)^r * choose(n, r)*(q/theta-r)^(n)/factorial(n)
return(sum(s))
}
})
pirwinhall_zero <- Vectorize(function(q, n, theta, p){
# function to calculate cumulative density of sum of n unif(0, theta) RVs
# inputs : q - quantile
# : n - number of unif RVs to sum
# : theta - population maximum
# outputs : numeric value that is the cumulative density at x
if(q/theta < 0) return(0)
if(q/theta > n) return(1)
if(q/theta >= 0 & q/theta <= n){
X <- floor(q/theta)
r <- seq(from = 0, to = X)
s <- (-1)^r * choose(n, r)*(q/theta-r)^(n)/factorial(n)
return(sum(s) - p)
}
})
qirwinhall <- Vectorize(function(p, n, theta){
# function to calculate quantile of sum of n unif(0, theta) RVs
# inputs : x - probability
# : n - number of unif RVs to sum
# : theta - population maximum
# outputs : numeric value that is the cumulative density at x
tmp <- uniroot(pirwinhall_zero, n = n, theta = theta, p = p, lower = 0, upper = n*theta, tol = .000000001)
return(tmp[[1]])
})
unif_sum_pwrfunc_greater <- function(theta, alpha, theta_not, n){
if(is.na(theta[length(theta)])){
return(NA)
}
k <- qirwinhall(p = 1 - alpha, n = n, theta = theta_not)
return(1 - pirwinhall(q = k, n = n, theta = theta))
}
unif_sum_pwrfunc_less <- function(theta, alpha, theta_not, n){
if(is.na(theta[length(theta)])){
return(NA)
}
k <- qirwinhall(p = alpha, n = n, theta = theta_not)
return(pirwinhall(q = k, n = n, theta = theta))
}
unif_sum_pwrfunc_noteqto <- function(theta, alpha, theta_not, n){
if(is.na(theta[length(theta)])){
return(NA)
}
k1 <- qirwinhall(p = alpha/2, n = n, theta = theta_not)
k2 <- qirwinhall(p = 1 - alpha/2, n = n, theta = theta_not)
return(1 - pirwinhall(q = k2, n = n, theta = theta) + pirwinhall(q = k1, n = n, theta = theta))
}
##############################
### SAMPLING DISTRIBUTIONS ###
##############################
# exponential
exp.samp <- function(statistic, alternative, theta, theta.not, n, alpha){
if(statistic == "Sum of the X's"){
if(alternative == 'Not equal to'){
plot(1, type = 'n', axes = F, xlab = '', ylab = '')
text(1, 1, 'Sampling distributions for not equal to alternatives are still under development.', col = 2)
} else if(alternative == 'Greater than'){
plot(1, type = "n", las = 1,
xlim = c(0.001, max(c(2*n*theta, 2*n*theta.not))),
ylim = c(0, max(
c(max(dgamma(seq(0, 2*n*theta), n, 1 / theta)),
max(dgamma(seq(0, 2*n*theta.not), n, 1 / theta.not)))
)),
xlab = "T(x)",
ylab = bquote(f[Sigma(X[i])](x)),
main = bquote("Sampling Distribution of T(X) ="~Sigma(X[i])~"for"~theta~"="~.(round(theta, 2))~"and"~theta[0]~"="~.(round(theta.not,2))))
(alt.curve <- curve(dgamma(x, n, 1 / theta), add = T, n = 1000))
(null.curve <- curve(dgamma(x, n, 1 / theta.not), add = T, lty = 2, n = 1000))
g.star <- qgamma(alpha, n, 1 / theta.not, lower.tail = F)
abline(v = g.star, lty = 4, col = "red")
if(1 - pchisq(theta.not*qchisq(1 - alpha, 2*n)/theta, 2*n) >= alpha ) {
polygon(x = c(g.star, g.star, alt.curve$x[which(alt.curve$x > g.star)], g.star),
y = c(0, dgamma(g.star, n, 1 / theta),
dgamma(alt.curve$x[which(alt.curve$x > g.star)], n, 1/theta), 0),
col = "grey")
polygon(x = c(g.star, g.star, null.curve$x[which(null.curve$x > g.star)], g.star),
y = c(0, dgamma(g.star, n, 1 / theta.not),
dgamma(null.curve$x[which(null.curve$x > g.star)], n, 1/theta.not), 0),
col = "red")
} else {
polygon(x = c(g.star, g.star, null.curve$x[which(null.curve$x > g.star)], g.star),
y = c(0, dgamma(g.star, n, 1 / theta.not),
dgamma(null.curve$x[which(null.curve$x > g.star)], n, 1/theta.not), 0),
col = "red")
polygon(x = c(g.star, g.star, alt.curve$x[which(alt.curve$x > g.star)], g.star),
y = c(0, dgamma(g.star, n, 1 / theta),
dgamma(alt.curve$x[which(alt.curve$x > g.star)], n, 1/theta), 0),
col = "grey")
}
power <- round(1 - pchisq(theta.not*qchisq(1 - alpha, 2*n)/theta, 2*n), 3)
legend("topright",
legend = c(expression(paste("Sampling Distribution Under ",theta)), bquote("Sampling Distribution Under"~theta[0]), bquote(alpha), bquote(beta(theta)~"="~.(power))),
lty = c(1,2,NA,NA), pch = c(NA, NA, 15, 15), col = c(1,1,"red","gray") , bty = "n")
text(x = g.star, y = dgamma(g.star, n, 1 / theta.not), labels = bquote('crit val'~'='~.(round(g.star, 2))), col = "red", pos = 2)
} else if(alternative == 'Less than'){
plot(1, type = "n", las = 1,
xlim = c(0.001, max(c(2*n*theta, 2*n*theta.not))),
ylim = c(0, max(
c(max(dgamma(seq(0, 2*n*theta), n, 1 / theta)),
max(dgamma(seq(0, 2*n*theta.not), n, 1 / theta.not)))
)),
xlab = "T(x)",
ylab = bquote(f[Sigma(X[i])](x)),
main = bquote("Sampling Distribution of T(X) ="~Sigma(X[i])~"for"~theta~"="~.(round(theta, 2))~"and"~theta[0]~"="~.(round(theta.not,2))))
(alt.curve <- curve(dgamma(x, n, 1 / theta), add = T, n = 1000))
(null.curve <- curve(dgamma(x, n, 1 / theta.not), add = T, lty = 2, n = 1000))
g.star <- qgamma(alpha, n, 1 / theta.not, lower.tail = T)
abline(v = g.star, lty = 4, col = "red")
if(pchisq(theta.not*qchisq(alpha, 2*n)/theta, 2*n) >= alpha ) {
polygon(x = c(g.star, alt.curve$x[which(alt.curve$x < g.star)], g.star, g.star),
y = c(0, dgamma(alt.curve$x[which(alt.curve$x < g.star)], n, 1/theta),
dgamma(g.star, n, 1 / theta), 0),
col = "grey")
polygon(x = c(g.star, null.curve$x[which(null.curve$x < g.star)], g.star, g.star),
y = c(0, dgamma(null.curve$x[which(null.curve$x < g.star)], n, 1/theta.not),
dgamma(g.star, n, 1 / theta.not), 0),
col = "red")
} else {
polygon(x = c(g.star, null.curve$x[which(null.curve$x < g.star)], g.star, g.star),
y = c(0, dgamma(null.curve$x[which(null.curve$x < g.star)], n, 1/theta.not),
dgamma(g.star, n, 1 / theta.not), 0),
col = "red")
polygon(x = c(g.star, alt.curve$x[which(alt.curve$x < g.star)], g.star, g.star),
y = c(0, dgamma(alt.curve$x[which(alt.curve$x < g.star)], n, 1/theta),
dgamma(g.star, n, 1 / theta), 0),
col = "grey")
}
legend("topright",
legend = c(expression(paste("Sampling Distribution Under ",theta)), bquote("Sampling Distribution Under"~theta[0]), bquote(alpha), bquote(beta(theta)~"="~.(round(pchisq(theta.not*qchisq(alpha, 2*n)/theta, 2*n), 3)))),
lty = c(1,2,NA,NA), pch = c(NA, NA, 15, 15), col = c(1,1,"red","gray") , bty = "n")
#text(x = 100, y = .001, labels = paste("Power = ", dgamma(theta, n, 1 / theta)))
text(x = g.star, y = dgamma(g.star, n, 1 / theta.not), labels = bquote('crit val'~'='~.(round(g.star, 2))), col = "red", pos = 4)
}
}
if(statistic == "Sample Minimum"){
if(alternative == "Not equal to"){
plot(1, type = 'n', axes = F, xlab = '', ylab = '')
text(1, 1, 'Sampling distributions for not equal to alternatives are still under development.', col = 2)
} else if(alternative == 'Less than'){
plot(1, type = "n", las = 1,
xlim = c(0.001, max(c(5*theta/n, 5*theta.not/n))), #5 times gives at least 96% of sampling dist by Chebychevs
ylim = c(0, max(
c(max(dexp(seq(0, 5*theta/n), n / theta)),
max(dexp(seq(0, 5*theta.not/n), n / theta.not)))
)),
xlab = "T(x)",
ylab = bquote(f[X[(1)]](x)),
main = bquote("Sampling Distribution of T(X) ="~X[(1)]~"for"~theta~"="~.(round(theta, 2))~"and"~theta[0]~"="~.(round(theta.not,2))))
(alt.curve <- curve(dexp(x, n / theta), add = T, n = 1000))
(null.curve <- curve(dexp(x, n / theta.not), add = T, lty = 2, n = 1000))
g.star <- qexp(alpha, n / theta.not, lower.tail = T)
abline(v = g.star, lty = 4, col = "red")
if(1 - exp(theta.not*log(1 - alpha)/theta) >= alpha ) {
polygon(x = c(0, 0, alt.curve$x[which(alt.curve$x < g.star)], g.star, g.star),
y = c(0, dexp(0, n/ theta), dexp(alt.curve$x[which(alt.curve$x < g.star)], n /theta),
dexp(g.star, n/ theta), 0), col = "grey")
polygon(x = c(0, null.curve$x[which(null.curve$x < g.star)], g.star, g.star),
y = c(0, dexp(null.curve$x[which(null.curve$x < g.star)], n / theta.not),
dexp(g.star, n / theta.not), 0), col = "red")
} else {
polygon(x = c(0, null.curve$x[which(null.curve$x < g.star)], g.star, g.star),
y = c(0, dexp(null.curve$x[which(null.curve$x < g.star)], n / theta.not),
dexp(g.star, n / theta.not), 0), col = "red")
polygon(x = c(0, 0, alt.curve$x[which(alt.curve$x < g.star)], g.star, g.star),
y = c(0, dexp(0, n/ theta), dexp(alt.curve$x[which(alt.curve$x < g.star)], n /theta),
dexp(g.star, n/ theta), 0), col = "grey")
}
legend("topright",
legend = c(expression(paste("Sampling Distribution Under ",theta)), bquote("Sampling Distribution Under"~theta[0]), bquote(alpha), bquote(beta(theta)~"="~.(round(1 - exp(theta.not*log(1 - alpha)/theta), 2)))),
lty = c(1,2,NA,NA), pch = c(NA, NA, 15, 15), col = c(1,1,"red","gray") , bty = "n")
text(x = g.star, y = dgamma(g.star, n, 1 / theta.not), labels = bquote('crit val'~'='~.(round(g.star, 2))), col = "red", pos = 4)
} else if(alternative == 'Greater than'){
plot(1, type = "n", las = 1,
xlim = c(0.001, max(c(5*theta/n, 5*theta.not/n))), #5 times gives at least 96.5% of sampling dist by Chebychevs
ylim = c(0, max(
c(max(dexp(seq(0, 5*theta/n), n / theta)),
max(dexp(seq(0, 5*theta.not/n), n / theta.not)))
)),
xlab = "T(x)",
ylab = bquote(f[X[(1)]](x)),
main = bquote("Sampling Distribution of T(X) ="~X[(1)]~"for"~theta~"="~.(round(theta, 2))~"and"~theta[0]~"="~.(round(theta.not,2))))
(alt.curve <- curve(dexp(x, n / theta), add = T, n = 1000))
(null.curve <- curve(dexp(x, n / theta.not), add = T, lty = 2, n = 1000))
g.star <- qexp(alpha, n / theta.not, lower.tail = F)
abline(v = g.star, lty = 4, col = "red")
if(exp(theta.not*log(alpha)/theta) >= alpha ) {
polygon(x = c(g.star, g.star, alt.curve$x[which(alt.curve$x > g.star)], g.star),
y = c(0, dexp(g.star, n/ theta),
dexp(alt.curve$x[which(alt.curve$x > g.star)], n /theta), 0),
col = "grey")
polygon(x = c(g.star, g.star, null.curve$x[which(null.curve$x > g.star)], g.star),
y = c(0, dexp(g.star, n / theta.not),
dexp(null.curve$x[which(null.curve$x > g.star)], n / theta.not), 0),
col = "red")
} else {
polygon(x = c(g.star, g.star, null.curve$x[which(null.curve$x > g.star)], g.star),
y = c(0, dexp(g.star, n / theta.not),
dexp(null.curve$x[which(null.curve$x > g.star)], n / theta.not), 0),
col = "red")
polygon(x = c(g.star, g.star, alt.curve$x[which(alt.curve$x > g.star)], g.star),
y = c(0, dexp(g.star, n/ theta),
dexp(alt.curve$x[which(alt.curve$x > g.star)], n /theta), 0),
col = "grey")
}
legend("topright",
legend = c(expression(paste("Sampling Distribution Under ",theta)), bquote("Sampling Distribution Under"~theta[0]), bquote(alpha), bquote(beta(theta)~"="~.(round(exp(theta.not*log(alpha)/theta), 2)))),
lty = c(1,2,NA,NA), pch = c(NA, NA, 15, 15), col = c(1,1,"red","gray") , bty = "n")
text(x = g.star, y = dexp(g.star, n / theta.not), labels = bquote('crit val'~'='~.(round(g.star, 2))), col = "red", pos = 2)
}
}
if(statistic == 'Sample Maximum'){
if(alternative == 'Not equal to'){
plot(1, type = 'n', axes = F, xlab = '', ylab = '')
text(1, 1, 'Sampling distributions for not equal to alternatives are still under development.', col = 2)
}
if(alternative == 'Less than'){
# sampling distribution functions
cdf <- function(x, y, theta, n){
(1 - exp(-x/theta))^n - y
}
sampdist <- function(x, theta, n){
n*(1 - exp(-x/theta))^(n-1) * exp(-x/theta) * 1/theta
}
# crit val
g.star <- -theta.not*log(1 - alpha^(1/n))
# plotting limits
upper.x <- max(c(
uniroot(cdf, y = .99, theta = theta, n = n, lower = 0, upper = 5*n*theta)$root,
uniroot(cdf, y = .99, theta = theta.not, n = n, lower = 0, upper = 5*n*theta.not)$root
))
nullmax <- max(sapply(seq(from = 0, to = upper.x, by = .01), FUN = function(x) sampdist(x, theta.not, n)))
altmax <- max(sapply(seq(from = 0, to = upper.x, by = .01), FUN = function(x) sampdist(x, theta, n)))
upper.y <- max(c(nullmax, altmax))
# plot
plot(1, type = "n", las = 1,
xlim = c(0.001, upper.x),
ylim = c(0, upper.y),
xlab = "T(x)",
ylab = bquote(f[X[(n)]](x)),
main = bquote("Sampling Distribution of T(X) ="~X[(n)]~"for"~theta~"="~.(round(theta, 2))~"and"~theta[0]~"="~.(round(theta.not,2))))
(alt.curve <- curve(sampdist(x, theta, n), add = T, n = 1000))
(null.curve <- curve(sampdist(x, theta.not, n), add = T, lty = 2, n = 1000))
abline(v = g.star, lty = 4, col = "red")
# polygons
if(exp_samp_max_pwrfunc_less(theta, alpha, theta.not, n) >= alpha){
polygon(x = c(0, 0, alt.curve$x[which(alt.curve$x < g.star)], g.star, g.star),
y = c(0, sampdist(0, theta, n), sampdist(alt.curve$x[which(alt.curve$x < g.star)], theta, n),
sampdist(g.star, theta, n), 0), col = "grey")
polygon(x = c(0, 0, alt.curve$x[which(alt.curve$x < g.star)], g.star, g.star),
y = c(0, sampdist(0, theta.not, n), sampdist(alt.curve$x[which(alt.curve$x < g.star)], theta.not, n),
sampdist(g.star, theta.not, n), 0), col = "red")
} else{
polygon(x = c(0, 0, alt.curve$x[which(alt.curve$x < g.star)], g.star, g.star),
y = c(0, sampdist(0, theta.not, n), sampdist(alt.curve$x[which(alt.curve$x < g.star)], theta.not, n),
sampdist(g.star, theta.not, n), 0), col = "red")
polygon(x = c(0, 0, alt.curve$x[which(alt.curve$x < g.star)], g.star, g.star),
y = c(0, sampdist(0, theta, n), sampdist(alt.curve$x[which(alt.curve$x < g.star)], theta, n),
sampdist(g.star, theta, n), 0), col = "grey")
}
# legend
power <- round(exp_samp_max_pwrfunc_less(theta, alpha, theta.not, n), 2)
legend("topright",
legend = c(expression(paste("Sampling Distribution Under ",theta)), bquote("Sampling Distribution Under"~theta[0]), bquote(alpha), bquote(beta(theta)~"="~.(power))),
lty = c(1,2,NA,NA), pch = c(NA, NA, 15, 15), col = c(1,1,"red","gray") , bty = "n")
text(x = g.star, y = sampdist(g.star, theta.not, n), labels = bquote('crit val'~'='~.(round(g.star, 2))), col = "red", pos = 4)
}
if(alternative == 'Greater than'){
# sampling distribution functions
cdf <- function(x, y, theta, n){
(1 - exp(-x/theta))^n - y
}
sampdist <- function(x, theta, n){
n*(1 - exp(-x/theta))^(n-1) * exp(-x/theta) * 1/theta
}
# crit val
g.star <- -theta.not*log(1 - (1 - alpha)^(1/n))
# plotting limits
upper.x <- max(c(
uniroot(cdf, y = .999, theta = theta, n = n, lower = 0, upper = 5*n*theta)$root,
uniroot(cdf, y = .999, theta = theta.not, n = n, lower = 0, upper = 5*n*theta.not)$root
))
nullmax <- max(sapply(seq(from = 0, to = upper.x, by = .01), FUN = function(x) sampdist(x, theta.not, n)))
altmax <- max(sapply(seq(from = 0, to = upper.x, by = .01), FUN = function(x) sampdist(x, theta, n)))
upper.y <- max(c(nullmax, altmax))
# plot
plot(1, type = "n", las = 1,
xlim = c(0.001, upper.x),
ylim = c(0, upper.y),
xlab = "T(x)",
ylab = bquote(f[X[(n)]](x)),
main = bquote("Sampling Distribution of T(X) ="~X[(n)]~"for"~theta~"="~.(round(theta, 2))~"and"~theta[0]~"="~.(round(theta.not,2))))
(alt.curve <- curve(sampdist(x, theta, n), add = T, n = 1000))
(null.curve <- curve(sampdist(x, theta.not, n), add = T, lty = 2, n = 1000))
abline(v = g.star, lty = 4, col = "red")
# polygons
if(exp_samp_max_pwrfunc_greater(theta, alpha, theta.not, n) >= alpha){
polygon(x = c(g.star, g.star, alt.curve$x[which(alt.curve$x > g.star)], g.star),
y = c(0, sampdist(g.star, theta, n), sampdist(alt.curve$x[which(alt.curve$x > g.star)], theta, n), 0),
col = "grey")
polygon(x = c(g.star, g.star, null.curve$x[which(null.curve$x > g.star)], g.star),
y = c(0, sampdist(g.star, theta.not, n),
sampdist(null.curve$x[which(null.curve$x > g.star)], theta.not, n), 0),
col = "red")
} else{
polygon(x = c(g.star, g.star, null.curve$x[which(null.curve$x > g.star)], g.star),
y = c(0, sampdist(g.star, theta.not, n),
sampdist(null.curve$x[which(null.curve$x > g.star)], theta.not, n), 0),
col = "red")
polygon(x = c(g.star, g.star, alt.curve$x[which(alt.curve$x > g.star)], g.star),
y = c(0, sampdist(g.star, theta, n), sampdist(alt.curve$x[which(alt.curve$x > g.star)], theta, n), 0),
col = "grey")
}
# legend
power <- round(exp_samp_max_pwrfunc_greater(theta, alpha, theta.not, n), 2)
legend("topright",
legend = c(expression(paste("Sampling Distribution Under ",theta)), bquote("Sampling Distribution Under"~theta[0]), bquote(alpha), bquote(beta(theta)~"="~.(power))),
lty = c(1,2,NA,NA), pch = c(NA, NA, 15, 15), col = c(1,1,"red","gray") , bty = "n")
text(x = g.star, y = sampdist(g.star, theta.not, n), labels = bquote('crit val'~'='~.(round(g.star, 2))), col = "red", pos = 2)
}
}
}
# normal
norm.samp <- function(statistic, alternative, theta, theta.not, n, alpha, sigma){
if(statistic == "Sum of the X's"){
if(alternative == 'Not equal to'){
plot(1, type = 'n', axes = F, xlab = '', ylab = '')
text(1, 1, 'Sampling distributions for not equal to alternatives are still under development.', col = 2)
}
if(alternative == 'Less than'){
# crit val
g.star <- qnorm(alpha, n*theta.not, sqrt(n)*sigma)
# plotting limits
upper.x <- max(c(
n*theta.not + 5*sqrt(n)*sigma,
n*theta + 5*sqrt(n)*sigma
))
lower.x <- min(c(
n*theta.not - 5*sqrt(n)*sigma,
n*theta - 5*sqrt(n)*sigma
))
nullmax <- max(sapply(seq(from = lower.x, to = upper.x, by = .01), FUN = function(x) dnorm(x, theta.not*n, sqrt(n)*sigma)))
altmax <- max(sapply(seq(from = lower.x, to = upper.x, by = .01), FUN = function(x) dnorm(x, theta*n, sqrt(n)*sigma)))
upper.y <- max(c(nullmax, altmax))
# plot
plot(1, type = "n", las = 1,
xlim = c(lower.x, upper.x),
ylim = c(0, upper.y),
xlab = "T(x)",
ylab = bquote(f[Sigma(X[i])](x)),
main = bquote("Sampling Distribution of T(X) ="~Sigma(X[i])~"for"~theta~"="~.(round(theta, 2))~"and"~theta[0]~"="~.(round(theta.not,2))))
(alt.curve <- curve(dnorm(x, n*theta, sqrt(n)*sigma), add = T, n = 1000))
(null.curve <- curve(dnorm(x, n*theta.not, sqrt(n)*sigma), add = T, lty = 2, n = 1000))
abline(v = g.star, lty = 4, col = "red")
# polygons
if(pnorm(qnorm(alpha, n*theta.not, sqrt(n)*sigma), n*theta, sqrt(n)*sigma) >= alpha){
polygon(x = c(alt.curve$x[which(alt.curve$x < g.star)], g.star, g.star),
y = c(dnorm(alt.curve$x[which(alt.curve$x < g.star)], n*theta, sqrt(n)*sigma),
dnorm(g.star, n*theta, sqrt(n)*sigma), 0), col = "grey")
polygon(x = c(alt.curve$x[which(alt.curve$x < g.star)], g.star, g.star),
y = c(dnorm(alt.curve$x[which(alt.curve$x < g.star)], n*theta.not, sqrt(n)*sigma),
dnorm(g.star, n*theta.not, sqrt(n)*sigma), 0), col = "red")
} else{
polygon(x = c(alt.curve$x[which(alt.curve$x < g.star)], g.star, g.star),
y = c(dnorm(alt.curve$x[which(alt.curve$x < g.star)], n*theta.not, sqrt(n)*sigma),
dnorm(g.star, n*theta.not, sqrt(n)*sigma), 0), col = "red")
polygon(x = c(alt.curve$x[which(alt.curve$x < g.star)], g.star, g.star),
y = c(dnorm(alt.curve$x[which(alt.curve$x < g.star)], n*theta, sqrt(n)*sigma),
dnorm(g.star, n*theta, sqrt(n)*sigma), 0), col = "grey")
}
# legend
power <- round(pnorm(qnorm(alpha, n*theta.not, sqrt(n)*sigma), n*theta, sqrt(n)*sigma), 2)
legend("topright",
legend = c(expression(paste("Sampling Distribution Under ",theta)), bquote("Sampling Distribution Under"~theta[0]), bquote(alpha), bquote(beta(theta)~"="~.(power))),
lty = c(1,2,NA,NA), pch = c(NA, NA, 15, 15), col = c(1,1,"red","gray") , bty = "n")
text(x = g.star, y = dnorm(g.star, n*theta.not, sqrt(n)*sigma), labels = bquote('crit val'~'='~.(round(g.star, 2))), col = "red", pos = 4)
}
if(alternative == 'Greater than'){
# crit val
g.star <- qnorm(1 - alpha, n*theta.not, sqrt(n)*sigma)
# plotting limits
upper.x <- max(c(
n*theta.not + 5*sqrt(n)*sigma,
n*theta + 5*sqrt(n)*sigma
))
lower.x <- min(c(
n*theta.not - 5*sqrt(n)*sigma,
n*theta - 5*sqrt(n)*sigma
))
nullmax <- max(sapply(seq(from = lower.x, to = upper.x, by = .01), FUN = function(x) dnorm(x, theta.not*n, sqrt(n)*sigma)))
altmax <- max(sapply(seq(from = lower.x, to = upper.x, by = .01), FUN = function(x) dnorm(x, theta*n, sqrt(n)*sigma)))
upper.y <- max(c(nullmax, altmax))
# plot
plot(1, type = "n", las = 1,
xlim = c(lower.x, upper.x),
ylim = c(0, upper.y),
xlab = "T(x)",
ylab = bquote(f[Sigma(X[i])](x)),
main = bquote("Sampling Distribution of T(X) ="~Sigma(X[i])~"for"~theta~"="~.(round(theta, 2))~"and"~theta[0]~"="~.(round(theta.not,2))))
(alt.curve <- curve(dnorm(x, n*theta, sqrt(n)*sigma), add = T, n = 1000))
(null.curve <- curve(dnorm(x, n*theta.not, sqrt(n)*sigma), add = T, lty = 2, n = 1000))
abline(v = g.star, lty = 4, col = "red")
# polygons
if(1 - pnorm(qnorm(1 - alpha, n*theta.not, sqrt(n)*sigma), n*theta, sqrt(n)*sigma) >= alpha){
polygon(x = c(g.star, g.star, alt.curve$x[which(alt.curve$x > g.star)]),
y = c(0, dnorm(g.star, theta*n, sqrt(n)*sigma), dnorm(alt.curve$x[which(alt.curve$x > g.star)], n*theta, sqrt(n)*sigma)),
col = "grey")
polygon(x = c(g.star, g.star, null.curve$x[which(null.curve$x > g.star)]),
y = c(0, dnorm(g.star, n*theta.not, sqrt(n)*sigma),
dnorm(null.curve$x[which(null.curve$x > g.star)], n*theta.not, sqrt(n)*sigma)),
col = "red")
} else{
polygon(x = c(g.star, g.star, null.curve$x[which(null.curve$x > g.star)]),
y = c(0, dnorm(g.star, n*theta.not, sqrt(n)*sigma),
dnorm(null.curve$x[which(null.curve$x > g.star)], n*theta.not, sqrt(n)*sigma)),
col = "red")
polygon(x = c(g.star, g.star, alt.curve$x[which(alt.curve$x > g.star)]),
y = c(0, dnorm(g.star, theta*n, sqrt(n)*sigma), dnorm(alt.curve$x[which(alt.curve$x > g.star)], n*theta, sqrt(n)*sigma)),
col = "grey")
}
# legend
power <- round(1 - pnorm(qnorm(1 - alpha, n*theta.not, sqrt(n)*sigma), n*theta, sqrt(n)*sigma), 2)
legend("topright",
legend = c(expression(paste("Sampling Distribution Under ",theta)), bquote("Sampling Distribution Under"~theta[0]), bquote(alpha), bquote(beta(theta)~"="~.(power))),
lty = c(1,2,NA,NA), pch = c(NA, NA, 15, 15), col = c(1,1,"red","gray") , bty = "n")
text(x = g.star, y = dnorm(g.star, n*theta.not, sqrt(n)*sigma), labels = bquote('crit val'~'='~.(round(g.star, 2))), col = "red", pos = 2)
}
}
if(statistic == "Sample Minimum"){
if(alternative == 'Not equal to'){
plot(1, type = 'n', axes = F, xlab = '', ylab = '')
text(1, 1, 'Sampling distributions for not equal to alternatives are still under development.', col = 2)
}
if(alternative == 'Less than'){
# sampling distribution functions
cdf <- function(x, y, theta, n, sigma = sigma){
1 - (1 - pnorm(x, theta, sigma))^n - y
}
sampdist <- function(x, theta, n, sigma = sigma){
n * (1 - pnorm(x, theta, sigma))^(n-1) * dnorm(x, theta, sigma)
}
# crit val
g.star <- qnorm(1 - (1 - alpha)^(1/n), theta.not, sigma)
# plotting limit
upper.x <- max(c(
uniroot(cdf, y = .9999, theta = theta, n = n, sigma = sigma, lower = theta - 5*n*sigma, upper = theta + 5*n*sigma)$root,
uniroot(cdf, y = .9999, theta = theta.not, n = n, sigma = sigma, lower = theta.not - 5*n*sigma, upper = theta.not + 5*n*sigma)$root
))
lower.x <- min(c(
uniroot(cdf, y = .0001, theta = theta, n = n, sigma = sigma, lower = theta - 5*n*sigma, upper = theta + 5*n*sigma)$root,
uniroot(cdf, y = .0001, theta = theta.not, n = n, sigma = sigma, lower = theta.not - 5*n*sigma, upper = theta.not + 5*n*sigma)$root
))
nullmax <- max(sapply(seq(from = lower.x, to = upper.x, by = .01), FUN = function(x) sampdist(x, theta.not, n, sigma)))
altmax <- max(sapply(seq(from = lower.x, to = upper.x, by = .01), FUN = function(x) sampdist(x, theta, n, sigma)))
upper.y <- max(c(nullmax, altmax))
# plot
plot(1, type = "n", las = 1,
xlim = c(lower.x, upper.x),
ylim = c(0, upper.y),
xlab = "T(x)",
ylab = bquote(f[X[(1)]](x)),
main = bquote("Sampling Distribution of T(X) ="~X[(1)]~"for"~theta~"="~.(round(theta, 2))~"and"~theta[0]~"="~.(round(theta.not,2))))
(alt.curve <- curve(sampdist(x, theta, n, sigma), add = T, n = 1000))
(null.curve <- curve(sampdist(x, theta.not, n, sigma), add = T, lty = 2, n = 1000))
abline(v = g.star, lty = 4, col = "red")
# polygons
if(norm_samp_min_pwrfunc_less(theta, alpha, sigma, theta.not, n) >= alpha){
polygon(x = c(alt.curve$x[which(alt.curve$x < g.star)], g.star, g.star),
y = c(sampdist(alt.curve$x[which(alt.curve$x < g.star)], theta, n, sigma),
sampdist(g.star, theta, n, sigma), 0), col = "grey")
polygon(x = c(alt.curve$x[which(alt.curve$x < g.star)], g.star, g.star),
y = c(sampdist(alt.curve$x[which(alt.curve$x < g.star)], theta.not, n, sigma),
sampdist(g.star, theta.not, n, sigma), 0), col = "red")
} else{
polygon(x = c(alt.curve$x[which(alt.curve$x < g.star)], g.star, g.star),
y = c(sampdist(alt.curve$x[which(alt.curve$x < g.star)], theta.not, n, sigma),
sampdist(g.star, theta.not, n, sigma), 0), col = "red")
polygon(x = c(alt.curve$x[which(alt.curve$x < g.star)], g.star, g.star),
y = c(sampdist(alt.curve$x[which(alt.curve$x < g.star)], theta, n, sigma),
sampdist(g.star, theta, n, sigma), 0), col = "grey")
}
# legend
power <- round(norm_samp_min_pwrfunc_less(theta, alpha, sigma, theta.not, n), 3)
legend("topright",
legend = c(expression(paste("Sampling Distribution Under ",theta)), bquote("Sampling Distribution Under"~theta[0]), bquote(alpha), bquote(beta(theta)~"="~.(power))),
lty = c(1,2,NA,NA), pch = c(NA, NA, 15, 15), col = c(1,1,"red","gray") , bty = "n")
text(x = g.star, y = sampdist(g.star, theta.not, n, sigma), labels = bquote('crit val'~'='~.(round(g.star, 2))), col = "red", pos = 4)
}
if(alternative == 'Greater than'){
# sampling distribution functions
cdf <- function(x, y, theta, n, sigma = sigma){
(pnorm(x, theta, sigma))^n - y
}
sampdist <- function(x, theta, n, sigma = sigma){
n* (pnorm(x, theta, sigma))^(n-1) * dnorm(x, theta, sigma)
}
# crit val
g.star <- qnorm((1 - alpha)^(1/n), theta.not, sigma)
# plotting limit
upper.x <- max(c(
uniroot(cdf, y = .9999, theta = theta, n = n, sigma = sigma, lower = theta - 5*n*sigma, upper = theta + 5*n*sigma)$root,
uniroot(cdf, y = .9999, theta = theta.not, n = n, sigma = sigma, lower = theta.not - 5*n*sigma, upper = theta.not + 5*n*sigma)$root
))
lower.x <- min(c(
uniroot(cdf, y = .0001, theta = theta, n = n, sigma = sigma, lower = theta - 5*n*sigma, upper = theta + 5*n*sigma)$root,
uniroot(cdf, y = .0001, theta = theta.not, n = n, sigma = sigma, lower = theta.not - 5*n*sigma, upper = theta.not + 5*n*sigma)$root
))
nullmax <- max(sapply(seq(from = lower.x, to = upper.x, by = .01), FUN = function(x) sampdist(x, theta.not, n, sigma)))
altmax <- max(sapply(seq(from = lower.x, to = upper.x, by = .01), FUN = function(x) sampdist(x, theta, n, sigma)))
upper.y <- max(c(nullmax, altmax))
# plot
plot(1, type = "n", las = 1,
xlim = c(lower.x, upper.x),
ylim = c(0, upper.y),
xlab = "T(x)",
ylab = bquote(f[X[(n)]](x)),
main = bquote("Sampling Distribution of T(X) ="~X[(n)]~"for"~theta~"="~.(round(theta, 2))~"and"~theta[0]~"="~.(round(theta.not,2))))
(alt.curve <- curve(sampdist(x, theta, n, sigma), add = T, n = 1000))
(null.curve <- curve(sampdist(x, theta.not, n, sigma), add = T, lty = 2, n = 1000))
abline(v = g.star, lty = 4, col = "red")
# polygons
if(norm_samp_max_pwrfunc_greater(theta, alpha, sigma, theta.not, n) >= alpha){
polygon(x = c(g.star, g.star, alt.curve$x[which(alt.curve$x > g.star)]),
y = c(0, sampdist(g.star, theta, n, sigma), sampdist(alt.curve$x[which(alt.curve$x > g.star)], theta, n, sigma)),
col = "grey")
polygon(x = c(g.star, g.star, null.curve$x[which(null.curve$x > g.star)]),
y = c(0, sampdist(g.star, theta.not, n, sigma),
sampdist(null.curve$x[which(null.curve$x > g.star)], theta.not, n, sigma)),
col = "red")
} else{
polygon(x = c(g.star, g.star, null.curve$x[which(null.curve$x > g.star)]),
y = c(0, sampdist(g.star, theta.not, n, sigma),
sampdist(null.curve$x[which(null.curve$x > g.star)], theta.not, n, sigma)),
col = "red")
polygon(x = c(g.star, g.star, alt.curve$x[which(alt.curve$x > g.star)]),
y = c(0, sampdist(g.star, theta, n, sigma), sampdist(alt.curve$x[which(alt.curve$x > g.star)], theta, n, sigma)),
col = "grey")
}
# legend
power <- round(norm_samp_min_pwrfunc_greater(theta, alpha, sigma, theta.not, n), 3)
legend("topright",
legend = c(expression(paste("Sampling Distribution Under ",theta)), bquote("Sampling Distribution Under"~theta[0]), bquote(alpha), bquote(beta(theta)~"="~.(power))),
lty = c(1,2,NA,NA), pch = c(NA, NA, 15, 15), col = c(1,1,"red","gray") , bty = "n")
text(x = g.star, y = sampdist(g.star, theta.not, n, sigma), labels = bquote('crit val'~'='~.(round(g.star, 2))), col = "red", pos = 2)
}
}
if(statistic == "Sample Maximum"){
if(alternative == 'Not equal to'){
plot(1, type = 'n', axes = F, xlab = '', ylab = '')
text(1, 1, 'Sampling distributions for not equal to alternatives are still under development.', col = 2)
}
if(alternative == 'Less than'){
# sampling distribution functions
cdf <- function(x, y, theta, n, sigma = sigma){
(pnorm(x, theta, sigma))^n - y
}
sampdist <- function(x, theta, n, sigma = sigma){
n* (pnorm(x, theta, sigma))^(n-1) * dnorm(x, theta, sigma)
}
# crit val
g.star <- qnorm(alpha^(1/n), theta.not, sigma)
# plotting limit
upper.x <- max(c(
uniroot(cdf, y = .9999, theta = theta, n = n, sigma = sigma, lower = theta - 5*n*sigma, upper = theta + 5*n*sigma)$root,
uniroot(cdf, y = .9999, theta = theta.not, n = n, sigma = sigma, lower = theta.not - 5*n*sigma, upper = theta.not + 5*n*sigma)$root
))
lower.x <- min(c(
uniroot(cdf, y = .0001, theta = theta, n = n, sigma = sigma, lower = theta - 5*n*sigma, upper = theta + 5*n*sigma)$root,
uniroot(cdf, y = .0001, theta = theta.not, n = n, sigma = sigma, lower = theta.not - 5*n*sigma, upper = theta.not + 5*n*sigma)$root
))
nullmax <- max(sapply(seq(from = lower.x, to = upper.x, by = .01), FUN = function(x) sampdist(x, theta.not, n, sigma)))
altmax <- max(sapply(seq(from = lower.x, to = upper.x, by = .01), FUN = function(x) sampdist(x, theta, n, sigma)))
upper.y <- max(c(nullmax, altmax))
# plot
plot(1, type = "n", las = 1,
xlim = c(lower.x, upper.x),
ylim = c(0, upper.y),
xlab = "T(x)",
ylab = bquote(f[X[(n)]](x)),
main = bquote("Sampling Distribution of T(X) ="~X[(n)]~"for"~theta~"="~.(round(theta, 2))~"and"~theta[0]~"="~.(round(theta.not,2))))
(alt.curve <- curve(sampdist(x, theta, n, sigma), add = T, n = 1000))
(null.curve <- curve(sampdist(x, theta.not, n, sigma), add = T, lty = 2, n = 1000))
abline(v = g.star, lty = 4, col = "red")
# polygons
if(norm_samp_max_pwrfunc_less(theta, alpha, sigma, theta.not, n) >= alpha){
polygon(x = c(alt.curve$x[which(alt.curve$x < g.star)], g.star, g.star),
y = c(sampdist(alt.curve$x[which(alt.curve$x < g.star)], theta, n, sigma),
sampdist(g.star, theta, n, sigma), 0), col = "grey")
polygon(x = c(alt.curve$x[which(alt.curve$x < g.star)], g.star, g.star),
y = c(sampdist(alt.curve$x[which(alt.curve$x < g.star)], theta.not, n, sigma),
sampdist(g.star, theta.not, n, sigma), 0), col = "red")
} else{
polygon(x = c(alt.curve$x[which(alt.curve$x < g.star)], g.star, g.star),
y = c(sampdist(alt.curve$x[which(alt.curve$x < g.star)], theta.not, n, sigma),
sampdist(g.star, theta.not, n, sigma), 0), col = "red")
polygon(x = c(alt.curve$x[which(alt.curve$x < g.star)], g.star, g.star),
y = c(sampdist(alt.curve$x[which(alt.curve$x < g.star)], theta, n, sigma),
sampdist(g.star, theta, n, sigma), 0), col = "grey")
}
# legend
power <- round(norm_samp_max_pwrfunc_less(theta, alpha, sigma, theta.not, n), 3)
legend("topright",
legend = c(expression(paste("Sampling Distribution Under ",theta)), bquote("Sampling Distribution Under"~theta[0]), bquote(alpha), bquote(beta(theta)~"="~.(power))),
lty = c(1,2,NA,NA), pch = c(NA, NA, 15, 15), col = c(1,1,"red","gray") , bty = "n")
text(x = g.star, y = sampdist(g.star, theta.not, n, sigma), labels = bquote('crit val'~'='~.(round(g.star, 2))), col = "red", pos = 4)
}
if(alternative == 'Greater than'){
# sampling distribution functions
cdf <- function(x, y, theta, n, sigma = sigma){
(pnorm(x, theta, sigma))^n - y
}
sampdist <- function(x, theta, n, sigma = sigma){
n* (pnorm(x, theta, sigma))^(n-1) * dnorm(x, theta, sigma)
}
# crit val
g.star <- qnorm((1 - alpha)^(1/n), theta.not, sigma)
# plotting limit
upper.x <- max(c(
uniroot(cdf, y = .9999, theta = theta, n = n, sigma = sigma, lower = theta - 5*n*sigma, upper = theta + 5*n*sigma)$root,
uniroot(cdf, y = .9999, theta = theta.not, n = n, sigma = sigma, lower = theta.not - 5*n*sigma, upper = theta.not + 5*n*sigma)$root
))
lower.x <- min(c(
uniroot(cdf, y = .0001, theta = theta, n = n, sigma = sigma, lower = theta - 5*n*sigma, upper = theta + 5*n*sigma)$root,
uniroot(cdf, y = .0001, theta = theta.not, n = n, sigma = sigma, lower = theta.not - 5*n*sigma, upper = theta.not + 5*n*sigma)$root
))
nullmax <- max(sapply(seq(from = lower.x, to = upper.x, by = .01), FUN = function(x) sampdist(x, theta.not, n, sigma)))
altmax <- max(sapply(seq(from = lower.x, to = upper.x, by = .01), FUN = function(x) sampdist(x, theta, n, sigma)))
upper.y <- max(c(nullmax, altmax))
# plot
plot(1, type = "n", las = 1,
xlim = c(lower.x, upper.x),
ylim = c(0, upper.y),
xlab = "T(x)",
ylab = bquote(f[X[(n)]](x)),
main = bquote("Sampling Distribution of T(X) ="~X[(n)]~"for"~theta~"="~.(round(theta, 2))~"and"~theta[0]~"="~.(round(theta.not,2))))
(alt.curve <- curve(sampdist(x, theta, n, sigma), add = T, n = 1000))
(null.curve <- curve(sampdist(x, theta.not, n, sigma), add = T, lty = 2, n = 1000))
abline(v = g.star, lty = 4, col = "red")
# polygons
if(norm_samp_max_pwrfunc_greater(theta, alpha, sigma, theta.not, n) >= alpha){
polygon(x = c(g.star, g.star, alt.curve$x[which(alt.curve$x > g.star)]),
y = c(0, sampdist(g.star, theta, n, sigma), sampdist(alt.curve$x[which(alt.curve$x > g.star)], theta, n, sigma)),
col = "grey")
polygon(x = c(g.star, g.star, null.curve$x[which(null.curve$x > g.star)]),
y = c(0, sampdist(g.star, theta.not, n, sigma),
sampdist(null.curve$x[which(null.curve$x > g.star)], theta.not, n, sigma)),
col = "red")
} else{
polygon(x = c(g.star, g.star, null.curve$x[which(null.curve$x > g.star)]),
y = c(0, sampdist(g.star, theta.not, n, sigma),
sampdist(null.curve$x[which(null.curve$x > g.star)], theta.not, n, sigma)),
col = "red")
polygon(x = c(g.star, g.star, alt.curve$x[which(alt.curve$x > g.star)]),
y = c(0, sampdist(g.star, theta, n, sigma), sampdist(alt.curve$x[which(alt.curve$x > g.star)], theta, n, sigma)),
col = "grey")
}
# legend
power <- round(norm_samp_max_pwrfunc_greater(theta, alpha, sigma, theta.not, n), 3)
legend("topright",
legend = c(expression(paste("Sampling Distribution Under ",theta)), bquote("Sampling Distribution Under"~theta[0]), bquote(alpha), bquote(beta(theta)~"="~.(power))),
lty = c(1,2,NA,NA), pch = c(NA, NA, 15, 15), col = c(1,1,"red","gray") , bty = "n")
text(x = g.star, y = sampdist(g.star, theta.not, n, sigma), labels = bquote('crit val'~'='~.(round(g.star, 2))), col = "red", pos = 2)
}
}
}
# exponential
unif.samp <- function(statistic, alternative, theta, theta.not, n, alpha){
if(statistic == "Sum of the X's"){
if(alternative == 'Not equal to'){
plot(1, type = 'n', axes = F, xlab = '', ylab = '')
text(1, 1, 'Sampling distributions for not equal to alternatives are still under development.', col = 2)
}
if(alternative == 'Less than'){
# crit val
g.star <- qirwinhall(alpha, n, theta.not)
sampdist <- function(x, theta, n, sigma = sigma){
n * (1 - pnorm(x, theta, sigma))^(n-1) * dnorm(x, theta, sigma)
}
# plotting limit
upper.x <- max(c(
uniroot(pirwinhall_zero, p = .99, theta = theta, n = n, lower = theta - 5*n, upper = theta + 5*n)$root,
uniroot(pirwinhall_zero, p = .99, theta = theta.not, n = n, lower = theta.not - 5*n, upper = theta.not + 5*n)$root
))
lower.x <- min(c(
uniroot(pirwinhall_zero, p = .01, theta = theta, n = n, lower = theta - 5*n, upper = theta + 5*n)$root,
uniroot(pirwinhall_zero, p = .01, theta = theta.not, n = n, lower = theta.not - 5*n, upper = theta.not + 5*n)$root
))
#################
### FIX THIS ####
#################
nullmax <- max(sapply(seq(from = lower.x, to = upper.x, by = .01), FUN = function(x) dnorm(x, theta.not*n, sqrt(n)*sigma)))
altmax <- max(sapply(seq(from = lower.x, to = upper.x, by = .01), FUN = function(x) dnorm(x, theta*n, sqrt(n)*sigma)))
upper.y <- max(c(nullmax, altmax))
# plot
plot(1, type = "n", las = 1,
xlim = c(lower.x, upper.x),
ylim = c(0, upper.y),
xlab = "T(x)",
ylab = bquote(f[Sigma(X[i])](x)),
main = bquote("Sampling Distribution of T(X) ="~Sigma(X[i])~"for"~theta~"="~.(round(theta, 2))~"and"~theta[0]~"="~.(round(theta.not,2))))
(alt.curve <- curve(dnorm(x, n*theta, sqrt(n)*sigma), add = T, n = 1000))
(null.curve <- curve(dnorm(x, n*theta.not, sqrt(n)*sigma), add = T, lty = 2, n = 1000))
abline(v = g.star, lty = 4, col = "red")
# polygons
if(pnorm(qnorm(alpha, n*theta.not, sqrt(n)*sigma), n*theta, sqrt(n)*sigma) >= alpha){
polygon(x = c(alt.curve$x[which(alt.curve$x < g.star)], g.star, g.star),
y = c(dnorm(alt.curve$x[which(alt.curve$x < g.star)], n*theta, sqrt(n)*sigma),
dnorm(g.star, n*theta, sqrt(n)*sigma), 0), col = "grey")
polygon(x = c(alt.curve$x[which(alt.curve$x < g.star)], g.star, g.star),
y = c(dnorm(alt.curve$x[which(alt.curve$x < g.star)], n*theta.not, sqrt(n)*sigma),
dnorm(g.star, n*theta.not, sqrt(n)*sigma), 0), col = "red")
} else{
polygon(x = c(alt.curve$x[which(alt.curve$x < g.star)], g.star, g.star),
y = c(dnorm(alt.curve$x[which(alt.curve$x < g.star)], n*theta.not, sqrt(n)*sigma),
dnorm(g.star, n*theta.not, sqrt(n)*sigma), 0), col = "red")
polygon(x = c(alt.curve$x[which(alt.curve$x < g.star)], g.star, g.star),
y = c(dnorm(alt.curve$x[which(alt.curve$x < g.star)], n*theta, sqrt(n)*sigma),
dnorm(g.star, n*theta, sqrt(n)*sigma), 0), col = "grey")
}
# legend
power <- round(pnorm(qnorm(alpha, n*theta.not, sqrt(n)*sigma), n*theta, sqrt(n)*sigma), 2)
legend("topright",
legend = c(expression(paste("Sampling Distribution Under ",theta)), bquote("Sampling Distribution Under"~theta[0]), bquote(alpha), bquote(beta(theta)~"="~.(power))),
lty = c(1,2,NA,NA), pch = c(NA, NA, 15, 15), col = c(1,1,"red","gray") , bty = "n")
text(x = g.star, y = dnorm(g.star, n*theta.not, sqrt(n)*sigma), labels = bquote('crit val'~'='~.(round(g.star, 2))), col = "red", pos = 4)
}
if(alternative == 'Greater than'){
}
}
if(statistic == "Sample Minimum"){
if(alternative == 'Not equal to'){
plot(1, type = 'n', axes = F, xlab = '', ylab = '')
text(1, 1, 'Sampling distributions for not equal to alternatives are still under development.', col = 2)
}
if(alternative == 'Less than'){
}
if(alternative == 'Greater than'){
}
}
if(statistic == "Sample Maximum"){
if(alternative == 'Not equal to'){
plot(1, type = 'n', axes = F, xlab = '', ylab = '')
text(1, 1, 'Sampling distributions for not equal to alternatives are still under development.', col = 2)
}
if(alternative == 'Less than'){
NULL
}
if(alternative == 'Greater than'){
NULL
}
}
}
StrattonCh/powerapppackage documentation built on May 21, 2019, 1:42 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
exp_samp_max_pwrfunc_greater <- function(theta, alpha, theta_not, n){
#function to calculate power function of sample max from exponential(theta) distribution
1 - (1 - exp(theta_not*log(1 - (1 - alpha)^(1/n)) / (theta) ))^n
}
exp_samp_max_pwrfunc_less <- function(theta, alpha, theta_not, n){
#function to calculate power function of sample max from exponential(theta) distribution
(1 - exp(theta_not*log(1 - (alpha)^(1/n)) / (theta) ))^n
}
exp_samp_max_pwrfunc_noteqto <- function(theta, alpha, theta_not, n){
#function to calculate power function of sample max from exponential(theta) distribution
1 - (1 - exp(theta_not*log(1 - (1 - alpha/2)^(1/n)) / (theta) ))^n + (1 - exp(theta_not*log(1 - (alpha/2)^(1/n)) / (theta) ))^n
}
norm_samp_min_pwrfunc_greater <- function(theta, alpha, sigma, theta_not, n){
#function to calculate power of sample min of normal(theta, sigma^2) distribution
k <- qnorm(1 - alpha^(1/n), theta_not, sigma)
(1 - pnorm(k, theta, sigma))^n
}
norm_samp_min_pwrfunc_less <- function(theta, alpha, sigma, theta_not, n){
#function to calculate power of sample min of normal(theta, sigma^2) distribution
k <- qnorm(1 - (1 - alpha)^(1/n), theta_not, sigma)
1 - (1 - pnorm(k, theta, sigma))^n
}
norm_min_cdf <- function(x, theta, sigma, n){
#function to calculate cdf of sample min from normal(theta, sigma^2)
1 - (1 - pnorm(x, theta, sigma))^n
}
norm_samp_min_pwrfunc_noteqto <- function(theta, alpha, sigma, theta_not, n){
#function to calculate power of sample min of normal(theta, sigma^2) distribution
k1 <- qnorm(1 - (1 - alpha/2)^(1/n), theta_not, sigma)
k2 <- qnorm(1 - (alpha/2)^(1/n), theta_not, sigma)
1 - norm_min_cdf(k2, theta, sigma, n) + norm_min_cdf(k1, theta, sigma, n)
}
norm_max_cdf <- function(x, theta, sigma, n){
#function to calculate the cdf of the sample max from normal(theta, sigma^2)
(pnorm(x, theta, sigma))^n
}
norm_samp_max_pwrfunc_greater <- function(theta, alpha, sigma, theta_not, n){
#function to calculate power of sample max of normal(theta, sigma^2) distribution
k <- qnorm((1 - alpha)^(1/n), theta_not, sigma)
1 - norm_max_cdf(k, theta, sigma, n)
}
norm_samp_max_pwrfunc_less <- function(theta, alpha, sigma, theta_not, n){
#function to calculate power of sample max of normal(theta, sigma^2) distribution
k <- qnorm(alpha^(1/n), theta_not, sigma)
norm_max_cdf(k, theta, sigma, n)
}
norm_samp_max_pwrfunc_noteqto <- function(theta, alpha, sigma, theta_not, n){
#function to calculate power of sample max of normal(theta, sigma^2) distribution
k1 <- qnorm((alpha/2)^(1/n), theta_not, sigma)
k2 <- qnorm((1 -alpha/2)^(1/n), theta_not, sigma)
1 - norm_max_cdf(k2, theta, sigma, n) + norm_max_cdf(k1, theta, sigma, n)
}
unif_samp_max_pwrfunc_greater <- Vectorize(function(theta, alpha, theta_not, n){
#function to calculate power of sample max of unif(0,theta)
k <- theta_not*(1 - alpha)^(1/n)
if(is.na(theta)){
return(NA)
}
if(theta < k){
return(0)
}
if(k <= theta){
return(1 - k^n/theta^n)
}
})
unif_samp_max_pwrfunc_less <- Vectorize(function(theta, alpha, theta_not, n){
#function to calculate power of sample max of unif(0,theta)
k <- theta_not*(alpha)^(1/n)
if(is.na(theta)){
return(NA)
}
if(theta < k){
return(1)
}
if(k <= theta){
return(k^n/theta^n)
}
})
unif_samp_max_pwrfunc_noteqto <- Vectorize(function(theta, alpha, theta_not, n){
#function to calculate power of sample max of unif(0,theta)
k1 <- theta_not*(alpha/2)^(1/n)
k2 <- theta_not*(1 - alpha/2)^(1/n)
if(is.na(theta)){
return(NA)
}
if(theta <= k1){
return(1)
}
if(k1 < theta & theta <= k2){
return(alpha*theta_not^n / (2*theta^n))
}
if(k2 <= theta){
return(1 - theta_not^n * (1 - alpha/2) / theta^n + alpha*theta_not^n / (2*theta^n))
}
})
unif_samp_min_pwrfunc_greater <- Vectorize(function(theta, alpha, theta_not, n){
#function to calculate power of sample min of unif(0,theta)
k <- theta_not*(1 - alpha^(1/n))
if(is.na(theta)){
return(NA)
}
if(theta <= k){
return(0)
}
if(theta > k){
return((1 - k/theta)^n)
}
})
unif_samp_min_pwrfunc_less <- Vectorize(function(theta, alpha, theta_not, n){
#function to calculate power of sample min of unif(0,theta)
k <- theta_not*(1 - (1 - alpha)^(1/n))
if(is.na(theta)){
return(NA)
}
if(theta <= k){
return(1)
}
if(theta > k){
return(1 - (1 - k/theta)^n)
}
})
unif_samp_min_pwrfunc_noteqto <- Vectorize(function(theta, alpha, theta_not, n){
#function to calculate power of sample min of unif(0,theta)
k1 <- theta_not*(1 - (1 - alpha/2)^(1/n))
k2 <- theta_not*(1 - (alpha/2)^(1/n))
if(is.na(theta)){
return(NA)
}
if(theta <= k1){
return(1)
}
if(k1 < theta & theta <= k2){
return(1 - (1 - k1/theta)^n)
}
if(theta > k2){
return(1 + (1 - k2/theta)^n - (1 - k1/theta)^n)
}
})
dirwinhall <- Vectorize(function(x, n, theta){
# function to calculate density of sum of n unif(0, theta) RVs
# inputs : x - value at which to calculate density
# : n - number of unif RVs to sum
# : theta - population maximum
# outputs : numeric value that is the density at x
if(x/theta < 0) return(0)
if(x/theta > n) return(0)
if(x/theta >= 0 & x/theta <= n){
X <- floor(x/theta)
r <- seq(from = 0, to = X)
s <- (-1)^r * choose(n, r)*(x/theta-r)^(n-1)/factorial(n-1)
return(sum(s)/theta)
}
})
pirwinhall <- Vectorize(function(q, n, theta){
# function to calculate cumulative density of sum of n unif(0, theta) RVs
# inputs : q - quantile
# : n - number of unif RVs to sum
# : theta - population maximum
# outputs : numeric value that is the cumulative density at x
if(q/theta < 0) return(0)
if(q/theta > n) return(1)
if(q/theta >= 0 & q/theta <= n){
X <- floor(q/theta)
r <- seq(from = 0, to = X)
s <- (-1)^r * choose(n, r)*(q/theta-r)^(n)/factorial(n)
return(sum(s))
}
})
pirwinhall_zero <- Vectorize(function(q, n, theta, p){
# function to calculate cumulative density of sum of n unif(0, theta) RVs
# inputs : q - quantile
# : n - number of unif RVs to sum
# : theta - population maximum
# outputs : numeric value that is the cumulative density at x
if(q/theta < 0) return(0)
if(q/theta > n) return(1)
if(q/theta >= 0 & q/theta <= n){
X <- floor(q/theta)
r <- seq(from = 0, to = X)
s <- (-1)^r * choose(n, r)*(q/theta-r)^(n)/factorial(n)
return(sum(s) - p)
}
})
qirwinhall <- Vectorize(function(p, n, theta){
# function to calculate quantile of sum of n unif(0, theta) RVs
# inputs : x - probability
# : n - number of unif RVs to sum
# : theta - population maximum
# outputs : numeric value that is the cumulative density at x
tmp <- uniroot(pirwinhall_zero, n = n, theta = theta, p = p, lower = 0, upper = n*theta, tol = .000000001)
return(tmp[[1]])
})
unif_sum_pwrfunc_greater <- function(theta, alpha, theta_not, n){
if(is.na(theta[length(theta)])){
return(NA)
}
k <- qirwinhall(p = 1 - alpha, n = n, theta = theta_not)
return(1 - pirwinhall(q = k, n = n, theta = theta))
}
unif_sum_pwrfunc_less <- function(theta, alpha, theta_not, n){
if(is.na(theta[length(theta)])){
return(NA)
}
k <- qirwinhall(p = alpha, n = n, theta = theta_not)
return(pirwinhall(q = k, n = n, theta = theta))
}
unif_sum_pwrfunc_noteqto <- function(theta, alpha, theta_not, n){
if(is.na(theta[length(theta)])){
return(NA)
}
k1 <- qirwinhall(p = alpha/2, n = n, theta = theta_not)
k2 <- qirwinhall(p = 1 - alpha/2, n = n, theta = theta_not)
return(1 - pirwinhall(q = k2, n = n, theta = theta) + pirwinhall(q = k1, n = n, theta = theta))
}
##############################
### SAMPLING DISTRIBUTIONS ###
##############################
# exponential
exp.samp <- function(statistic, alternative, theta, theta.not, n, alpha){
if(statistic == "Sum of the X's"){
if(alternative == 'Not equal to'){
plot(1, type = 'n', axes = F, xlab = '', ylab = '')
text(1, 1, 'Sampling distributions for not equal to alternatives are still under development.', col = 2)
} else if(alternative == 'Greater than'){
plot(1, type = "n", las = 1,
xlim = c(0.001, max(c(2*n*theta, 2*n*theta.not))),
ylim = c(0, max(
c(max(dgamma(seq(0, 2*n*theta), n, 1 / theta)),
max(dgamma(seq(0, 2*n*theta.not), n, 1 / theta.not)))
)),
xlab = "T(x)",
ylab = bquote(f[Sigma(X[i])](x)),
main = bquote("Sampling Distribution of T(X) ="~Sigma(X[i])~"for"~theta~"="~.(round(theta, 2))~"and"~theta[0]~"="~.(round(theta.not,2))))
(alt.curve <- curve(dgamma(x, n, 1 / theta), add = T, n = 1000))
(null.curve <- curve(dgamma(x, n, 1 / theta.not), add = T, lty = 2, n = 1000))
g.star <- qgamma(alpha, n, 1 / theta.not, lower.tail = F)
abline(v = g.star, lty = 4, col = "red")
if(1 - pchisq(theta.not*qchisq(1 - alpha, 2*n)/theta, 2*n) >= alpha ) {
polygon(x = c(g.star, g.star, alt.curve$x[which(alt.curve$x > g.star)], g.star),
y = c(0, dgamma(g.star, n, 1 / theta),
dgamma(alt.curve$x[which(alt.curve$x > g.star)], n, 1/theta), 0),
col = "grey")
polygon(x = c(g.star, g.star, null.curve$x[which(null.curve$x > g.star)], g.star),
y = c(0, dgamma(g.star, n, 1 / theta.not),
dgamma(null.curve$x[which(null.curve$x > g.star)], n, 1/theta.not), 0),
col = "red")
} else {
polygon(x = c(g.star, g.star, null.curve$x[which(null.curve$x > g.star)], g.star),
y = c(0, dgamma(g.star, n, 1 / theta.not),
dgamma(null.curve$x[which(null.curve$x > g.star)], n, 1/theta.not), 0),
col = "red")
polygon(x = c(g.star, g.star, alt.curve$x[which(alt.curve$x > g.star)], g.star),
y = c(0, dgamma(g.star, n, 1 / theta),
dgamma(alt.curve$x[which(alt.curve$x > g.star)], n, 1/theta), 0),
col = "grey")
}
power <- round(1 - pchisq(theta.not*qchisq(1 - alpha, 2*n)/theta, 2*n), 3)
legend("topright",
legend = c(expression(paste("Sampling Distribution Under ",theta)), bquote("Sampling Distribution Under"~theta[0]), bquote(alpha), bquote(beta(theta)~"="~.(power))),
lty = c(1,2,NA,NA), pch = c(NA, NA, 15, 15), col = c(1,1,"red","gray") , bty = "n")
text(x = g.star, y = dgamma(g.star, n, 1 / theta.not), labels = bquote('crit val'~'='~.(round(g.star, 2))), col = "red", pos = 2)
} else if(alternative == 'Less than'){
plot(1, type = "n", las = 1,
xlim = c(0.001, max(c(2*n*theta, 2*n*theta.not))),
ylim = c(0, max(
c(max(dgamma(seq(0, 2*n*theta), n, 1 / theta)),
max(dgamma(seq(0, 2*n*theta.not), n, 1 / theta.not)))
)),
xlab = "T(x)",
ylab = bquote(f[Sigma(X[i])](x)),
main = bquote("Sampling Distribution of T(X) ="~Sigma(X[i])~"for"~theta~"="~.(round(theta, 2))~"and"~theta[0]~"="~.(round(theta.not,2))))
(alt.curve <- curve(dgamma(x, n, 1 / theta), add = T, n = 1000))
(null.curve <- curve(dgamma(x, n, 1 / theta.not), add = T, lty = 2, n = 1000))
g.star <- qgamma(alpha, n, 1 / theta.not, lower.tail = T)
abline(v = g.star, lty = 4, col = "red")
if(pchisq(theta.not*qchisq(alpha, 2*n)/theta, 2*n) >= alpha ) {
polygon(x = c(g.star, alt.curve$x[which(alt.curve$x < g.star)], g.star, g.star),
y = c(0, dgamma(alt.curve$x[which(alt.curve$x < g.star)], n, 1/theta),
dgamma(g.star, n, 1 / theta), 0),
col = "grey")
polygon(x = c(g.star, null.curve$x[which(null.curve$x < g.star)], g.star, g.star),
y = c(0, dgamma(null.curve$x[which(null.curve$x < g.star)], n, 1/theta.not),
dgamma(g.star, n, 1 / theta.not), 0),
col = "red")
} else {
polygon(x = c(g.star, null.curve$x[which(null.curve$x < g.star)], g.star, g.star),
y = c(0, dgamma(null.curve$x[which(null.curve$x < g.star)], n, 1/theta.not),
dgamma(g.star, n, 1 / theta.not), 0),
col = "red")
polygon(x = c(g.star, alt.curve$x[which(alt.curve$x < g.star)], g.star, g.star),
y = c(0, dgamma(alt.curve$x[which(alt.curve$x < g.star)], n, 1/theta),
dgamma(g.star, n, 1 / theta), 0),
col = "grey")
}
legend("topright",
legend = c(expression(paste("Sampling Distribution Under ",theta)), bquote("Sampling Distribution Under"~theta[0]), bquote(alpha), bquote(beta(theta)~"="~.(round(pchisq(theta.not*qchisq(alpha, 2*n)/theta, 2*n), 3)))),
lty = c(1,2,NA,NA), pch = c(NA, NA, 15, 15), col = c(1,1,"red","gray") , bty = "n")
#text(x = 100, y = .001, labels = paste("Power = ", dgamma(theta, n, 1 / theta)))
text(x = g.star, y = dgamma(g.star, n, 1 / theta.not), labels = bquote('crit val'~'='~.(round(g.star, 2))), col = "red", pos = 4)
}
}
if(statistic == "Sample Minimum"){
if(alternative == "Not equal to"){
plot(1, type = 'n', axes = F, xlab = '', ylab = '')
text(1, 1, 'Sampling distributions for not equal to alternatives are still under development.', col = 2)
} else if(alternative == 'Less than'){
plot(1, type = "n", las = 1,
xlim = c(0.001, max(c(5*theta/n, 5*theta.not/n))), #5 times gives at least 96% of sampling dist by Chebychevs
ylim = c(0, max(
c(max(dexp(seq(0, 5*theta/n), n / theta)),
max(dexp(seq(0, 5*theta.not/n), n / theta.not)))
)),
xlab = "T(x)",
ylab = bquote(f[X[(1)]](x)),
main = bquote("Sampling Distribution of T(X) ="~X[(1)]~"for"~theta~"="~.(round(theta, 2))~"and"~theta[0]~"="~.(round(theta.not,2))))
(alt.curve <- curve(dexp(x, n / theta), add = T, n = 1000))
(null.curve <- curve(dexp(x, n / theta.not), add = T, lty = 2, n = 1000))
g.star <- qexp(alpha, n / theta.not, lower.tail = T)
abline(v = g.star, lty = 4, col = "red")
if(1 - exp(theta.not*log(1 - alpha)/theta) >= alpha ) {
polygon(x = c(0, 0, alt.curve$x[which(alt.curve$x < g.star)], g.star, g.star),
y = c(0, dexp(0, n/ theta), dexp(alt.curve$x[which(alt.curve$x < g.star)], n /theta),
dexp(g.star, n/ theta), 0), col = "grey")
polygon(x = c(0, null.curve$x[which(null.curve$x < g.star)], g.star, g.star),
y = c(0, dexp(null.curve$x[which(null.curve$x < g.star)], n / theta.not),
dexp(g.star, n / theta.not), 0), col = "red")
} else {
polygon(x = c(0, null.curve$x[which(null.curve$x < g.star)], g.star, g.star),
y = c(0, dexp(null.curve$x[which(null.curve$x < g.star)], n / theta.not),
dexp(g.star, n / theta.not), 0), col = "red")
polygon(x = c(0, 0, alt.curve$x[which(alt.curve$x < g.star)], g.star, g.star),
y = c(0, dexp(0, n/ theta), dexp(alt.curve$x[which(alt.curve$x < g.star)], n /theta),
dexp(g.star, n/ theta), 0), col = "grey")
}
legend("topright",
legend = c(expression(paste("Sampling Distribution Under ",theta)), bquote("Sampling Distribution Under"~theta[0]), bquote(alpha), bquote(beta(theta)~"="~.(round(1 - exp(theta.not*log(1 - alpha)/theta), 2)))),
lty = c(1,2,NA,NA), pch = c(NA, NA, 15, 15), col = c(1,1,"red","gray") , bty = "n")
text(x = g.star, y = dgamma(g.star, n, 1 / theta.not), labels = bquote('crit val'~'='~.(round(g.star, 2))), col = "red", pos = 4)
} else if(alternative == 'Greater than'){
plot(1, type = "n", las = 1,
xlim = c(0.001, max(c(5*theta/n, 5*theta.not/n))), #5 times gives at least 96.5% of sampling dist by Chebychevs
ylim = c(0, max(
c(max(dexp(seq(0, 5*theta/n), n / theta)),
max(dexp(seq(0, 5*theta.not/n), n / theta.not)))
)),
xlab = "T(x)",
ylab = bquote(f[X[(1)]](x)),
main = bquote("Sampling Distribution of T(X) ="~X[(1)]~"for"~theta~"="~.(round(theta, 2))~"and"~theta[0]~"="~.(round(theta.not,2))))
(alt.curve <- curve(dexp(x, n / theta), add = T, n = 1000))
(null.curve <- curve(dexp(x, n / theta.not), add = T, lty = 2, n = 1000))
g.star <- qexp(alpha, n / theta.not, lower.tail = F)
abline(v = g.star, lty = 4, col = "red")
if(exp(theta.not*log(alpha)/theta) >= alpha ) {
polygon(x = c(g.star, g.star, alt.curve$x[which(alt.curve$x > g.star)], g.star),
y = c(0, dexp(g.star, n/ theta),
dexp(alt.curve$x[which(alt.curve$x > g.star)], n /theta), 0),
col = "grey")
polygon(x = c(g.star, g.star, null.curve$x[which(null.curve$x > g.star)], g.star),
y = c(0, dexp(g.star, n / theta.not),
dexp(null.curve$x[which(null.curve$x > g.star)], n / theta.not), 0),
col = "red")
} else {
polygon(x = c(g.star, g.star, null.curve$x[which(null.curve$x > g.star)], g.star),
y = c(0, dexp(g.star, n / theta.not),
dexp(null.curve$x[which(null.curve$x > g.star)], n / theta.not), 0),
col = "red")
polygon(x = c(g.star, g.star, alt.curve$x[which(alt.curve$x > g.star)], g.star),
y = c(0, dexp(g.star, n/ theta),
dexp(alt.curve$x[which(alt.curve$x > g.star)], n /theta), 0),
col = "grey")
}
legend("topright",
legend = c(expression(paste("Sampling Distribution Under ",theta)), bquote("Sampling Distribution Under"~theta[0]), bquote(alpha), bquote(beta(theta)~"="~.(round(exp(theta.not*log(alpha)/theta), 2)))),
lty = c(1,2,NA,NA), pch = c(NA, NA, 15, 15), col = c(1,1,"red","gray") , bty = "n")
text(x = g.star, y = dexp(g.star, n / theta.not), labels = bquote('crit val'~'='~.(round(g.star, 2))), col = "red", pos = 2)
}
}
if(statistic == 'Sample Maximum'){
if(alternative == 'Not equal to'){
plot(1, type = 'n', axes = F, xlab = '', ylab = '')
text(1, 1, 'Sampling distributions for not equal to alternatives are still under development.', col = 2)
}
if(alternative == 'Less than'){
# sampling distribution functions
cdf <- function(x, y, theta, n){
(1 - exp(-x/theta))^n - y
}
sampdist <- function(x, theta, n){
n*(1 - exp(-x/theta))^(n-1) * exp(-x/theta) * 1/theta
}
# crit val
g.star <- -theta.not*log(1 - alpha^(1/n))
# plotting limits
upper.x <- max(c(
uniroot(cdf, y = .99, theta = theta, n = n, lower = 0, upper = 5*n*theta)$root,
uniroot(cdf, y = .99, theta = theta.not, n = n, lower = 0, upper = 5*n*theta.not)$root
))
nullmax <- max(sapply(seq(from = 0, to = upper.x, by = .01), FUN = function(x) sampdist(x, theta.not, n)))
altmax <- max(sapply(seq(from = 0, to = upper.x, by = .01), FUN = function(x) sampdist(x, theta, n)))
upper.y <- max(c(nullmax, altmax))
# plot
plot(1, type = "n", las = 1,
xlim = c(0.001, upper.x),
ylim = c(0, upper.y),
xlab = "T(x)",
ylab = bquote(f[X[(n)]](x)),
main = bquote("Sampling Distribution of T(X) ="~X[(n)]~"for"~theta~"="~.(round(theta, 2))~"and"~theta[0]~"="~.(round(theta.not,2))))
(alt.curve <- curve(sampdist(x, theta, n), add = T, n = 1000))
(null.curve <- curve(sampdist(x, theta.not, n), add = T, lty = 2, n = 1000))
abline(v = g.star, lty = 4, col = "red")
# polygons
if(exp_samp_max_pwrfunc_less(theta, alpha, theta.not, n) >= alpha){
polygon(x = c(0, 0, alt.curve$x[which(alt.curve$x < g.star)], g.star, g.star),
y = c(0, sampdist(0, theta, n), sampdist(alt.curve$x[which(alt.curve$x < g.star)], theta, n),
sampdist(g.star, theta, n), 0), col = "grey")
polygon(x = c(0, 0, alt.curve$x[which(alt.curve$x < g.star)], g.star, g.star),
y = c(0, sampdist(0, theta.not, n), sampdist(alt.curve$x[which(alt.curve$x < g.star)], theta.not, n),
sampdist(g.star, theta.not, n), 0), col = "red")
} else{
polygon(x = c(0, 0, alt.curve$x[which(alt.curve$x < g.star)], g.star, g.star),
y = c(0, sampdist(0, theta.not, n), sampdist(alt.curve$x[which(alt.curve$x < g.star)], theta.not, n),
sampdist(g.star, theta.not, n), 0), col = "red")
polygon(x = c(0, 0, alt.curve$x[which(alt.curve$x < g.star)], g.star, g.star),
y = c(0, sampdist(0, theta, n), sampdist(alt.curve$x[which(alt.curve$x < g.star)], theta, n),
sampdist(g.star, theta, n), 0), col = "grey")
}
# legend
power <- round(exp_samp_max_pwrfunc_less(theta, alpha, theta.not, n), 2)
legend("topright",
legend = c(expression(paste("Sampling Distribution Under ",theta)), bquote("Sampling Distribution Under"~theta[0]), bquote(alpha), bquote(beta(theta)~"="~.(power))),
lty = c(1,2,NA,NA), pch = c(NA, NA, 15, 15), col = c(1,1,"red","gray") , bty = "n")
text(x = g.star, y = sampdist(g.star, theta.not, n), labels = bquote('crit val'~'='~.(round(g.star, 2))), col = "red", pos = 4)
}
if(alternative == 'Greater than'){
# sampling distribution functions
cdf <- function(x, y, theta, n){
(1 - exp(-x/theta))^n - y
}
sampdist <- function(x, theta, n){
n*(1 - exp(-x/theta))^(n-1) * exp(-x/theta) * 1/theta
}
# crit val
g.star <- -theta.not*log(1 - (1 - alpha)^(1/n))
# plotting limits
upper.x <- max(c(
uniroot(cdf, y = .999, theta = theta, n = n, lower = 0, upper = 5*n*theta)$root,
uniroot(cdf, y = .999, theta = theta.not, n = n, lower = 0, upper = 5*n*theta.not)$root
))
nullmax <- max(sapply(seq(from = 0, to = upper.x, by = .01), FUN = function(x) sampdist(x, theta.not, n)))
altmax <- max(sapply(seq(from = 0, to = upper.x, by = .01), FUN = function(x) sampdist(x, theta, n)))
upper.y <- max(c(nullmax, altmax))
# plot
plot(1, type = "n", las = 1,
xlim = c(0.001, upper.x),
ylim = c(0, upper.y),
xlab = "T(x)",
ylab = bquote(f[X[(n)]](x)),
main = bquote("Sampling Distribution of T(X) ="~X[(n)]~"for"~theta~"="~.(round(theta, 2))~"and"~theta[0]~"="~.(round(theta.not,2))))
(alt.curve <- curve(sampdist(x, theta, n), add = T, n = 1000))
(null.curve <- curve(sampdist(x, theta.not, n), add = T, lty = 2, n = 1000))
abline(v = g.star, lty = 4, col = "red")
# polygons
if(exp_samp_max_pwrfunc_greater(theta, alpha, theta.not, n) >= alpha){
polygon(x = c(g.star, g.star, alt.curve$x[which(alt.curve$x > g.star)], g.star),
y = c(0, sampdist(g.star, theta, n), sampdist(alt.curve$x[which(alt.curve$x > g.star)], theta, n), 0),
col = "grey")
polygon(x = c(g.star, g.star, null.curve$x[which(null.curve$x > g.star)], g.star),
y = c(0, sampdist(g.star, theta.not, n),
sampdist(null.curve$x[which(null.curve$x > g.star)], theta.not, n), 0),
col = "red")
} else{
polygon(x = c(g.star, g.star, null.curve$x[which(null.curve$x > g.star)], g.star),
y = c(0, sampdist(g.star, theta.not, n),
sampdist(null.curve$x[which(null.curve$x > g.star)], theta.not, n), 0),
col = "red")
polygon(x = c(g.star, g.star, alt.curve$x[which(alt.curve$x > g.star)], g.star),
y = c(0, sampdist(g.star, theta, n), sampdist(alt.curve$x[which(alt.curve$x > g.star)], theta, n), 0),
col = "grey")
}
# legend
power <- round(exp_samp_max_pwrfunc_greater(theta, alpha, theta.not, n), 2)
legend("topright",
legend = c(expression(paste("Sampling Distribution Under ",theta)), bquote("Sampling Distribution Under"~theta[0]), bquote(alpha), bquote(beta(theta)~"="~.(power))),
lty = c(1,2,NA,NA), pch = c(NA, NA, 15, 15), col = c(1,1,"red","gray") , bty = "n")
text(x = g.star, y = sampdist(g.star, theta.not, n), labels = bquote('crit val'~'='~.(round(g.star, 2))), col = "red", pos = 2)
}
}
}
# normal
norm.samp <- function(statistic, alternative, theta, theta.not, n, alpha, sigma){
if(statistic == "Sum of the X's"){
if(alternative == 'Not equal to'){
plot(1, type = 'n', axes = F, xlab = '', ylab = '')
text(1, 1, 'Sampling distributions for not equal to alternatives are still under development.', col = 2)
}
if(alternative == 'Less than'){
# crit val
g.star <- qnorm(alpha, n*theta.not, sqrt(n)*sigma)
# plotting limits
upper.x <- max(c(
n*theta.not + 5*sqrt(n)*sigma,
n*theta + 5*sqrt(n)*sigma
))
lower.x <- min(c(
n*theta.not - 5*sqrt(n)*sigma,
n*theta - 5*sqrt(n)*sigma
))
nullmax <- max(sapply(seq(from = lower.x, to = upper.x, by = .01), FUN = function(x) dnorm(x, theta.not*n, sqrt(n)*sigma)))
altmax <- max(sapply(seq(from = lower.x, to = upper.x, by = .01), FUN = function(x) dnorm(x, theta*n, sqrt(n)*sigma)))
upper.y <- max(c(nullmax, altmax))
# plot
plot(1, type = "n", las = 1,
xlim = c(lower.x, upper.x),
ylim = c(0, upper.y),
xlab = "T(x)",
ylab = bquote(f[Sigma(X[i])](x)),
main = bquote("Sampling Distribution of T(X) ="~Sigma(X[i])~"for"~theta~"="~.(round(theta, 2))~"and"~theta[0]~"="~.(round(theta.not,2))))
(alt.curve <- curve(dnorm(x, n*theta, sqrt(n)*sigma), add = T, n = 1000))
(null.curve <- curve(dnorm(x, n*theta.not, sqrt(n)*sigma), add = T, lty = 2, n = 1000))
abline(v = g.star, lty = 4, col = "red")
# polygons
if(pnorm(qnorm(alpha, n*theta.not, sqrt(n)*sigma), n*theta, sqrt(n)*sigma) >= alpha){
polygon(x = c(alt.curve$x[which(alt.curve$x < g.star)], g.star, g.star),
y = c(dnorm(alt.curve$x[which(alt.curve$x < g.star)], n*theta, sqrt(n)*sigma),
dnorm(g.star, n*theta, sqrt(n)*sigma), 0), col = "grey")
polygon(x = c(alt.curve$x[which(alt.curve$x < g.star)], g.star, g.star),
y = c(dnorm(alt.curve$x[which(alt.curve$x < g.star)], n*theta.not, sqrt(n)*sigma),
dnorm(g.star, n*theta.not, sqrt(n)*sigma), 0), col = "red")
} else{
polygon(x = c(alt.curve$x[which(alt.curve$x < g.star)], g.star, g.star),
y = c(dnorm(alt.curve$x[which(alt.curve$x < g.star)], n*theta.not, sqrt(n)*sigma),
dnorm(g.star, n*theta.not, sqrt(n)*sigma), 0), col = "red")
polygon(x = c(alt.curve$x[which(alt.curve$x < g.star)], g.star, g.star),
y = c(dnorm(alt.curve$x[which(alt.curve$x < g.star)], n*theta, sqrt(n)*sigma),
dnorm(g.star, n*theta, sqrt(n)*sigma), 0), col = "grey")
}
# legend
power <- round(pnorm(qnorm(alpha, n*theta.not, sqrt(n)*sigma), n*theta, sqrt(n)*sigma), 2)
legend("topright",
legend = c(expression(paste("Sampling Distribution Under ",theta)), bquote("Sampling Distribution Under"~theta[0]), bquote(alpha), bquote(beta(theta)~"="~.(power))),
lty = c(1,2,NA,NA), pch = c(NA, NA, 15, 15), col = c(1,1,"red","gray") , bty = "n")
text(x = g.star, y = dnorm(g.star, n*theta.not, sqrt(n)*sigma), labels = bquote('crit val'~'='~.(round(g.star, 2))), col = "red", pos = 4)
}
if(alternative == 'Greater than'){
# crit val
g.star <- qnorm(1 - alpha, n*theta.not, sqrt(n)*sigma)
# plotting limits
upper.x <- max(c(
n*theta.not + 5*sqrt(n)*sigma,
n*theta + 5*sqrt(n)*sigma
))
lower.x <- min(c(
n*theta.not - 5*sqrt(n)*sigma,
n*theta - 5*sqrt(n)*sigma
))
nullmax <- max(sapply(seq(from = lower.x, to = upper.x, by = .01), FUN = function(x) dnorm(x, theta.not*n, sqrt(n)*sigma)))
altmax <- max(sapply(seq(from = lower.x, to = upper.x, by = .01), FUN = function(x) dnorm(x, theta*n, sqrt(n)*sigma)))
upper.y <- max(c(nullmax, altmax))
# plot
plot(1, type = "n", las = 1,
xlim = c(lower.x, upper.x),
ylim = c(0, upper.y),
xlab = "T(x)",
ylab = bquote(f[Sigma(X[i])](x)),
main = bquote("Sampling Distribution of T(X) ="~Sigma(X[i])~"for"~theta~"="~.(round(theta, 2))~"and"~theta[0]~"="~.(round(theta.not,2))))
(alt.curve <- curve(dnorm(x, n*theta, sqrt(n)*sigma), add = T, n = 1000))
(null.curve <- curve(dnorm(x, n*theta.not, sqrt(n)*sigma), add = T, lty = 2, n = 1000))
abline(v = g.star, lty = 4, col = "red")
# polygons
if(1 - pnorm(qnorm(1 - alpha, n*theta.not, sqrt(n)*sigma), n*theta, sqrt(n)*sigma) >= alpha){
polygon(x = c(g.star, g.star, alt.curve$x[which(alt.curve$x > g.star)]),
y = c(0, dnorm(g.star, theta*n, sqrt(n)*sigma), dnorm(alt.curve$x[which(alt.curve$x > g.star)], n*theta, sqrt(n)*sigma)),
col = "grey")
polygon(x = c(g.star, g.star, null.curve$x[which(null.curve$x > g.star)]),
y = c(0, dnorm(g.star, n*theta.not, sqrt(n)*sigma),
dnorm(null.curve$x[which(null.curve$x > g.star)], n*theta.not, sqrt(n)*sigma)),
col = "red")
} else{
polygon(x = c(g.star, g.star, null.curve$x[which(null.curve$x > g.star)]),
y = c(0, dnorm(g.star, n*theta.not, sqrt(n)*sigma),
dnorm(null.curve$x[which(null.curve$x > g.star)], n*theta.not, sqrt(n)*sigma)),
col = "red")
polygon(x = c(g.star, g.star, alt.curve$x[which(alt.curve$x > g.star)]),
y = c(0, dnorm(g.star, theta*n, sqrt(n)*sigma), dnorm(alt.curve$x[which(alt.curve$x > g.star)], n*theta, sqrt(n)*sigma)),
col = "grey")
}
# legend
power <- round(1 - pnorm(qnorm(1 - alpha, n*theta.not, sqrt(n)*sigma), n*theta, sqrt(n)*sigma), 2)
legend("topright",
legend = c(expression(paste("Sampling Distribution Under ",theta)), bquote("Sampling Distribution Under"~theta[0]), bquote(alpha), bquote(beta(theta)~"="~.(power))),
lty = c(1,2,NA,NA), pch = c(NA, NA, 15, 15), col = c(1,1,"red","gray") , bty = "n")
text(x = g.star, y = dnorm(g.star, n*theta.not, sqrt(n)*sigma), labels = bquote('crit val'~'='~.(round(g.star, 2))), col = "red", pos = 2)
}
}
if(statistic == "Sample Minimum"){
if(alternative == 'Not equal to'){
plot(1, type = 'n', axes = F, xlab = '', ylab = '')
text(1, 1, 'Sampling distributions for not equal to alternatives are still under development.', col = 2)
}
if(alternative == 'Less than'){
# sampling distribution functions
cdf <- function(x, y, theta, n, sigma = sigma){
1 - (1 - pnorm(x, theta, sigma))^n - y
}
sampdist <- function(x, theta, n, sigma = sigma){
n * (1 - pnorm(x, theta, sigma))^(n-1) * dnorm(x, theta, sigma)
}
# crit val
g.star <- qnorm(1 - (1 - alpha)^(1/n), theta.not, sigma)
# plotting limit
upper.x <- max(c(
uniroot(cdf, y = .9999, theta = theta, n = n, sigma = sigma, lower = theta - 5*n*sigma, upper = theta + 5*n*sigma)$root,
uniroot(cdf, y = .9999, theta = theta.not, n = n, sigma = sigma, lower = theta.not - 5*n*sigma, upper = theta.not + 5*n*sigma)$root
))
lower.x <- min(c(
uniroot(cdf, y = .0001, theta = theta, n = n, sigma = sigma, lower = theta - 5*n*sigma, upper = theta + 5*n*sigma)$root,
uniroot(cdf, y = .0001, theta = theta.not, n = n, sigma = sigma, lower = theta.not - 5*n*sigma, upper = theta.not + 5*n*sigma)$root
))
nullmax <- max(sapply(seq(from = lower.x, to = upper.x, by = .01), FUN = function(x) sampdist(x, theta.not, n, sigma)))
altmax <- max(sapply(seq(from = lower.x, to = upper.x, by = .01), FUN = function(x) sampdist(x, theta, n, sigma)))
upper.y <- max(c(nullmax, altmax))
# plot
plot(1, type = "n", las = 1,
xlim = c(lower.x, upper.x),
ylim = c(0, upper.y),
xlab = "T(x)",
ylab = bquote(f[X[(1)]](x)),
main = bquote("Sampling Distribution of T(X) ="~X[(1)]~"for"~theta~"="~.(round(theta, 2))~"and"~theta[0]~"="~.(round(theta.not,2))))
(alt.curve <- curve(sampdist(x, theta, n, sigma), add = T, n = 1000))
(null.curve <- curve(sampdist(x, theta.not, n, sigma), add = T, lty = 2, n = 1000))
abline(v = g.star, lty = 4, col = "red")
# polygons
if(norm_samp_min_pwrfunc_less(theta, alpha, sigma, theta.not, n) >= alpha){
polygon(x = c(alt.curve$x[which(alt.curve$x < g.star)], g.star, g.star),
y = c(sampdist(alt.curve$x[which(alt.curve$x < g.star)], theta, n, sigma),
sampdist(g.star, theta, n, sigma), 0), col = "grey")
polygon(x = c(alt.curve$x[which(alt.curve$x < g.star)], g.star, g.star),
y = c(sampdist(alt.curve$x[which(alt.curve$x < g.star)], theta.not, n, sigma),
sampdist(g.star, theta.not, n, sigma), 0), col = "red")
} else{
polygon(x = c(alt.curve$x[which(alt.curve$x < g.star)], g.star, g.star),
y = c(sampdist(alt.curve$x[which(alt.curve$x < g.star)], theta.not, n, sigma),
sampdist(g.star, theta.not, n, sigma), 0), col = "red")
polygon(x = c(alt.curve$x[which(alt.curve$x < g.star)], g.star, g.star),
y = c(sampdist(alt.curve$x[which(alt.curve$x < g.star)], theta, n, sigma),
sampdist(g.star, theta, n, sigma), 0), col = "grey")
}
# legend
power <- round(norm_samp_min_pwrfunc_less(theta, alpha, sigma, theta.not, n), 3)
legend("topright",
legend = c(expression(paste("Sampling Distribution Under ",theta)), bquote("Sampling Distribution Under"~theta[0]), bquote(alpha), bquote(beta(theta)~"="~.(power))),
lty = c(1,2,NA,NA), pch = c(NA, NA, 15, 15), col = c(1,1,"red","gray") , bty = "n")
text(x = g.star, y = sampdist(g.star, theta.not, n, sigma), labels = bquote('crit val'~'='~.(round(g.star, 2))), col = "red", pos = 4)
}
if(alternative == 'Greater than'){
# sampling distribution functions
cdf <- function(x, y, theta, n, sigma = sigma){
(pnorm(x, theta, sigma))^n - y
}
sampdist <- function(x, theta, n, sigma = sigma){
n* (pnorm(x, theta, sigma))^(n-1) * dnorm(x, theta, sigma)
}
# crit val
g.star <- qnorm((1 - alpha)^(1/n), theta.not, sigma)
# plotting limit
upper.x <- max(c(
uniroot(cdf, y = .9999, theta = theta, n = n, sigma = sigma, lower = theta - 5*n*sigma, upper = theta + 5*n*sigma)$root,
uniroot(cdf, y = .9999, theta = theta.not, n = n, sigma = sigma, lower = theta.not - 5*n*sigma, upper = theta.not + 5*n*sigma)$root
))
lower.x <- min(c(
uniroot(cdf, y = .0001, theta = theta, n = n, sigma = sigma, lower = theta - 5*n*sigma, upper = theta + 5*n*sigma)$root,
uniroot(cdf, y = .0001, theta = theta.not, n = n, sigma = sigma, lower = theta.not - 5*n*sigma, upper = theta.not + 5*n*sigma)$root
))
nullmax <- max(sapply(seq(from = lower.x, to = upper.x, by = .01), FUN = function(x) sampdist(x, theta.not, n, sigma)))
altmax <- max(sapply(seq(from = lower.x, to = upper.x, by = .01), FUN = function(x) sampdist(x, theta, n, sigma)))
upper.y <- max(c(nullmax, altmax))
# plot
plot(1, type = "n", las = 1,
xlim = c(lower.x, upper.x),
ylim = c(0, upper.y),
xlab = "T(x)",
ylab = bquote(f[X[(n)]](x)),
main = bquote("Sampling Distribution of T(X) ="~X[(n)]~"for"~theta~"="~.(round(theta, 2))~"and"~theta[0]~"="~.(round(theta.not,2))))
(alt.curve <- curve(sampdist(x, theta, n, sigma), add = T, n = 1000))
(null.curve <- curve(sampdist(x, theta.not, n, sigma), add = T, lty = 2, n = 1000))
abline(v = g.star, lty = 4, col = "red")
# polygons
if(norm_samp_max_pwrfunc_greater(theta, alpha, sigma, theta.not, n) >= alpha){
polygon(x = c(g.star, g.star, alt.curve$x[which(alt.curve$x > g.star)]),
y = c(0, sampdist(g.star, theta, n, sigma), sampdist(alt.curve$x[which(alt.curve$x > g.star)], theta, n, sigma)),
col = "grey")
polygon(x = c(g.star, g.star, null.curve$x[which(null.curve$x > g.star)]),
y = c(0, sampdist(g.star, theta.not, n, sigma),
sampdist(null.curve$x[which(null.curve$x > g.star)], theta.not, n, sigma)),
col = "red")
} else{
polygon(x = c(g.star, g.star, null.curve$x[which(null.curve$x > g.star)]),
y = c(0, sampdist(g.star, theta.not, n, sigma),
sampdist(null.curve$x[which(null.curve$x > g.star)], theta.not, n, sigma)),
col = "red")
polygon(x = c(g.star, g.star, alt.curve$x[which(alt.curve$x > g.star)]),
y = c(0, sampdist(g.star, theta, n, sigma), sampdist(alt.curve$x[which(alt.curve$x > g.star)], theta, n, sigma)),
col = "grey")
}
# legend
power <- round(norm_samp_min_pwrfunc_greater(theta, alpha, sigma, theta.not, n), 3)
legend("topright",
legend = c(expression(paste("Sampling Distribution Under ",theta)), bquote("Sampling Distribution Under"~theta[0]), bquote(alpha), bquote(beta(theta)~"="~.(power))),
lty = c(1,2,NA,NA), pch = c(NA, NA, 15, 15), col = c(1,1,"red","gray") , bty = "n")
text(x = g.star, y = sampdist(g.star, theta.not, n, sigma), labels = bquote('crit val'~'='~.(round(g.star, 2))), col = "red", pos = 2)
}
}
if(statistic == "Sample Maximum"){
if(alternative == 'Not equal to'){
plot(1, type = 'n', axes = F, xlab = '', ylab = '')
text(1, 1, 'Sampling distributions for not equal to alternatives are still under development.', col = 2)
}
if(alternative == 'Less than'){
# sampling distribution functions
cdf <- function(x, y, theta, n, sigma = sigma){
(pnorm(x, theta, sigma))^n - y
}
sampdist <- function(x, theta, n, sigma = sigma){
n* (pnorm(x, theta, sigma))^(n-1) * dnorm(x, theta, sigma)
}
# crit val
g.star <- qnorm(alpha^(1/n), theta.not, sigma)
# plotting limit
upper.x <- max(c(
uniroot(cdf, y = .9999, theta = theta, n = n, sigma = sigma, lower = theta - 5*n*sigma, upper = theta + 5*n*sigma)$root,
uniroot(cdf, y = .9999, theta = theta.not, n = n, sigma = sigma, lower = theta.not - 5*n*sigma, upper = theta.not + 5*n*sigma)$root
))
lower.x <- min(c(
uniroot(cdf, y = .0001, theta = theta, n = n, sigma = sigma, lower = theta - 5*n*sigma, upper = theta + 5*n*sigma)$root,
uniroot(cdf, y = .0001, theta = theta.not, n = n, sigma = sigma, lower = theta.not - 5*n*sigma, upper = theta.not + 5*n*sigma)$root
))
nullmax <- max(sapply(seq(from = lower.x, to = upper.x, by = .01), FUN = function(x) sampdist(x, theta.not, n, sigma)))
altmax <- max(sapply(seq(from = lower.x, to = upper.x, by = .01), FUN = function(x) sampdist(x, theta, n, sigma)))
upper.y <- max(c(nullmax, altmax))
# plot
plot(1, type = "n", las = 1,
xlim = c(lower.x, upper.x),
ylim = c(0, upper.y),
xlab = "T(x)",
ylab = bquote(f[X[(n)]](x)),
main = bquote("Sampling Distribution of T(X) ="~X[(n)]~"for"~theta~"="~.(round(theta, 2))~"and"~theta[0]~"="~.(round(theta.not,2))))
(alt.curve <- curve(sampdist(x, theta, n, sigma), add = T, n = 1000))
(null.curve <- curve(sampdist(x, theta.not, n, sigma), add = T, lty = 2, n = 1000))
abline(v = g.star, lty = 4, col = "red")
# polygons
if(norm_samp_max_pwrfunc_less(theta, alpha, sigma, theta.not, n) >= alpha){
polygon(x = c(alt.curve$x[which(alt.curve$x < g.star)], g.star, g.star),
y = c(sampdist(alt.curve$x[which(alt.curve$x < g.star)], theta, n, sigma),
sampdist(g.star, theta, n, sigma), 0), col = "grey")
polygon(x = c(alt.curve$x[which(alt.curve$x < g.star)], g.star, g.star),
y = c(sampdist(alt.curve$x[which(alt.curve$x < g.star)], theta.not, n, sigma),
sampdist(g.star, theta.not, n, sigma), 0), col = "red")
} else{
polygon(x = c(alt.curve$x[which(alt.curve$x < g.star)], g.star, g.star),
y = c(sampdist(alt.curve$x[which(alt.curve$x < g.star)], theta.not, n, sigma),
sampdist(g.star, theta.not, n, sigma), 0), col = "red")
polygon(x = c(alt.curve$x[which(alt.curve$x < g.star)], g.star, g.star),
y = c(sampdist(alt.curve$x[which(alt.curve$x < g.star)], theta, n, sigma),
sampdist(g.star, theta, n, sigma), 0), col = "grey")
}
# legend
power <- round(norm_samp_max_pwrfunc_less(theta, alpha, sigma, theta.not, n), 3)
legend("topright",
legend = c(expression(paste("Sampling Distribution Under ",theta)), bquote("Sampling Distribution Under"~theta[0]), bquote(alpha), bquote(beta(theta)~"="~.(power))),
lty = c(1,2,NA,NA), pch = c(NA, NA, 15, 15), col = c(1,1,"red","gray") , bty = "n")
text(x = g.star, y = sampdist(g.star, theta.not, n, sigma), labels = bquote('crit val'~'='~.(round(g.star, 2))), col = "red", pos = 4)
}
if(alternative == 'Greater than'){
# sampling distribution functions
cdf <- function(x, y, theta, n, sigma = sigma){
(pnorm(x, theta, sigma))^n - y
}
sampdist <- function(x, theta, n, sigma = sigma){
n* (pnorm(x, theta, sigma))^(n-1) * dnorm(x, theta, sigma)
}
# crit val
g.star <- qnorm((1 - alpha)^(1/n), theta.not, sigma)
# plotting limit
upper.x <- max(c(
uniroot(cdf, y = .9999, theta = theta, n = n, sigma = sigma, lower = theta - 5*n*sigma, upper = theta + 5*n*sigma)$root,
uniroot(cdf, y = .9999, theta = theta.not, n = n, sigma = sigma, lower = theta.not - 5*n*sigma, upper = theta.not + 5*n*sigma)$root
))
lower.x <- min(c(
uniroot(cdf, y = .0001, theta = theta, n = n, sigma = sigma, lower = theta - 5*n*sigma, upper = theta + 5*n*sigma)$root,
uniroot(cdf, y = .0001, theta = theta.not, n = n, sigma = sigma, lower = theta.not - 5*n*sigma, upper = theta.not + 5*n*sigma)$root
))
nullmax <- max(sapply(seq(from = lower.x, to = upper.x, by = .01), FUN = function(x) sampdist(x, theta.not, n, sigma)))
altmax <- max(sapply(seq(from = lower.x, to = upper.x, by = .01), FUN = function(x) sampdist(x, theta, n, sigma)))
upper.y <- max(c(nullmax, altmax))
# plot
plot(1, type = "n", las = 1,
xlim = c(lower.x, upper.x),
ylim = c(0, upper.y),
xlab = "T(x)",
ylab = bquote(f[X[(n)]](x)),
main = bquote("Sampling Distribution of T(X) ="~X[(n)]~"for"~theta~"="~.(round(theta, 2))~"and"~theta[0]~"="~.(round(theta.not,2))))
(alt.curve <- curve(sampdist(x, theta, n, sigma), add = T, n = 1000))
(null.curve <- curve(sampdist(x, theta.not, n, sigma), add = T, lty = 2, n = 1000))
abline(v = g.star, lty = 4, col = "red")
# polygons
if(norm_samp_max_pwrfunc_greater(theta, alpha, sigma, theta.not, n) >= alpha){
polygon(x = c(g.star, g.star, alt.curve$x[which(alt.curve$x > g.star)]),
y = c(0, sampdist(g.star, theta, n, sigma), sampdist(alt.curve$x[which(alt.curve$x > g.star)], theta, n, sigma)),
col = "grey")
polygon(x = c(g.star, g.star, null.curve$x[which(null.curve$x > g.star)]),
y = c(0, sampdist(g.star, theta.not, n, sigma),
sampdist(null.curve$x[which(null.curve$x > g.star)], theta.not, n, sigma)),
col = "red")
} else{
polygon(x = c(g.star, g.star, null.curve$x[which(null.curve$x > g.star)]),
y = c(0, sampdist(g.star, theta.not, n, sigma),
sampdist(null.curve$x[which(null.curve$x > g.star)], theta.not, n, sigma)),
col = "red")
polygon(x = c(g.star, g.star, alt.curve$x[which(alt.curve$x > g.star)]),
y = c(0, sampdist(g.star, theta, n, sigma), sampdist(alt.curve$x[which(alt.curve$x > g.star)], theta, n, sigma)),
col = "grey")
}
# legend
power <- round(norm_samp_max_pwrfunc_greater(theta, alpha, sigma, theta.not, n), 3)
legend("topright",
legend = c(expression(paste("Sampling Distribution Under ",theta)), bquote("Sampling Distribution Under"~theta[0]), bquote(alpha), bquote(beta(theta)~"="~.(power))),
lty = c(1,2,NA,NA), pch = c(NA, NA, 15, 15), col = c(1,1,"red","gray") , bty = "n")
text(x = g.star, y = sampdist(g.star, theta.not, n, sigma), labels = bquote('crit val'~'='~.(round(g.star, 2))), col = "red", pos = 2)
}
}
}
# exponential
unif.samp <- function(statistic, alternative, theta, theta.not, n, alpha){
if(statistic == "Sum of the X's"){
if(alternative == 'Not equal to'){
plot(1, type = 'n', axes = F, xlab = '', ylab = '')
text(1, 1, 'Sampling distributions for not equal to alternatives are still under development.', col = 2)
}
if(alternative == 'Less than'){
# crit val
g.star <- qirwinhall(alpha, n, theta.not)
sampdist <- function(x, theta, n, sigma = sigma){
n * (1 - pnorm(x, theta, sigma))^(n-1) * dnorm(x, theta, sigma)
}
# plotting limit
upper.x <- max(c(
uniroot(pirwinhall_zero, p = .99, theta = theta, n = n, lower = theta - 5*n, upper = theta + 5*n)$root,
uniroot(pirwinhall_zero, p = .99, theta = theta.not, n = n, lower = theta.not - 5*n, upper = theta.not + 5*n)$root
))
lower.x <- min(c(
uniroot(pirwinhall_zero, p = .01, theta = theta, n = n, lower = theta - 5*n, upper = theta + 5*n)$root,
uniroot(pirwinhall_zero, p = .01, theta = theta.not, n = n, lower = theta.not - 5*n, upper = theta.not + 5*n)$root
))
#################
### FIX THIS ####
#################
nullmax <- max(sapply(seq(from = lower.x, to = upper.x, by = .01), FUN = function(x) dnorm(x, theta.not*n, sqrt(n)*sigma)))
altmax <- max(sapply(seq(from = lower.x, to = upper.x, by = .01), FUN = function(x) dnorm(x, theta*n, sqrt(n)*sigma)))
upper.y <- max(c(nullmax, altmax))
# plot
plot(1, type = "n", las = 1,
xlim = c(lower.x, upper.x),
ylim = c(0, upper.y),
xlab = "T(x)",
ylab = bquote(f[Sigma(X[i])](x)),
main = bquote("Sampling Distribution of T(X) ="~Sigma(X[i])~"for"~theta~"="~.(round(theta, 2))~"and"~theta[0]~"="~.(round(theta.not,2))))
(alt.curve <- curve(dnorm(x, n*theta, sqrt(n)*sigma), add = T, n = 1000))
(null.curve <- curve(dnorm(x, n*theta.not, sqrt(n)*sigma), add = T, lty = 2, n = 1000))
abline(v = g.star, lty = 4, col = "red")
# polygons
if(pnorm(qnorm(alpha, n*theta.not, sqrt(n)*sigma), n*theta, sqrt(n)*sigma) >= alpha){
polygon(x = c(alt.curve$x[which(alt.curve$x < g.star)], g.star, g.star),
y = c(dnorm(alt.curve$x[which(alt.curve$x < g.star)], n*theta, sqrt(n)*sigma),
dnorm(g.star, n*theta, sqrt(n)*sigma), 0), col = "grey")
polygon(x = c(alt.curve$x[which(alt.curve$x < g.star)], g.star, g.star),
y = c(dnorm(alt.curve$x[which(alt.curve$x < g.star)], n*theta.not, sqrt(n)*sigma),
dnorm(g.star, n*theta.not, sqrt(n)*sigma), 0), col = "red")
} else{
polygon(x = c(alt.curve$x[which(alt.curve$x < g.star)], g.star, g.star),
y = c(dnorm(alt.curve$x[which(alt.curve$x < g.star)], n*theta.not, sqrt(n)*sigma),
dnorm(g.star, n*theta.not, sqrt(n)*sigma), 0), col = "red")
polygon(x = c(alt.curve$x[which(alt.curve$x < g.star)], g.star, g.star),
y = c(dnorm(alt.curve$x[which(alt.curve$x < g.star)], n*theta, sqrt(n)*sigma),
dnorm(g.star, n*theta, sqrt(n)*sigma), 0), col = "grey")
}
# legend
power <- round(pnorm(qnorm(alpha, n*theta.not, sqrt(n)*sigma), n*theta, sqrt(n)*sigma), 2)
legend("topright",
legend = c(expression(paste("Sampling Distribution Under ",theta)), bquote("Sampling Distribution Under"~theta[0]), bquote(alpha), bquote(beta(theta)~"="~.(power))),
lty = c(1,2,NA,NA), pch = c(NA, NA, 15, 15), col = c(1,1,"red","gray") , bty = "n")
text(x = g.star, y = dnorm(g.star, n*theta.not, sqrt(n)*sigma), labels = bquote('crit val'~'='~.(round(g.star, 2))), col = "red", pos = 4)
}
if(alternative == 'Greater than'){
}
}
if(statistic == "Sample Minimum"){
if(alternative == 'Not equal to'){
plot(1, type = 'n', axes = F, xlab = '', ylab = '')
text(1, 1, 'Sampling distributions for not equal to alternatives are still under development.', col = 2)
}
if(alternative == 'Less than'){
}
if(alternative == 'Greater than'){
}
}
if(statistic == "Sample Maximum"){
if(alternative == 'Not equal to'){
plot(1, type = 'n', axes = F, xlab = '', ylab = '')
text(1, 1, 'Sampling distributions for not equal to alternatives are still under development.', col = 2)
}
if(alternative == 'Less than'){
NULL
}
if(alternative == 'Greater than'){
NULL
}
}
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.