library(ggplot2)
# standard normal distribution
x <- normalGrid(r = 3)
plot(x$z, x$wgts)
# verify that numerical integration replicates sigma
# get grid points and weights
x <- normalGrid(mu = 2, sigma = 3)
# compute squared deviation from mean for grid points
dev <- (x$z - 2)^2
# multiply squared deviations by integration weights and sum
sigma2 <- sum(dev * x$wgts)
# square root of sigma2 should be sigma (3)
sqrt(sigma2)
# do it again with larger r to increase accuracy
x <- normalGrid(r = 22, mu = 2, sigma = 3)
sqrt(sum((x$z - 2)^2 * x$wgts))
# this can also be done by combining gridwgts and density
sqrt(sum((x$z - 2)^2 * x$gridwgts * x$density))
# integrate normal density and compare to built-in function
# to compute probability of being within 1 standard deviation
# of the mean
pnorm(1) - pnorm(-1)
x <- normalGrid(bounds = c(-1, 1))
sum(x$wgts)
sum(x$gridwgts * x$density)
# find expected sample size for default design with
# n.fix=1000
x <- gsDesign(n.fix = 1000)
x
# set a prior distribution for theta
y <- normalGrid(r = 3, mu = x$theta[2], sigma = x$theta[2] / 1.5)
z <- gsProbability(
k = 3, theta = y$z, n.I = x$n.I, a = x$lower$bound,
b = x$upper$bound
)
z <- gsProbability(d = x, theta = y$z)
cat(
"Expected sample size averaged over normal\n prior distribution for theta with \n mu=",
x$theta[2], "sigma=", x$theta[2] / 1.5, ":",
round(sum(z$en * y$wgt), 1), "\n"
)
plot(y$z, z$en,
xlab = "theta", ylab = "E{N}",
main = "Expected sample size for different theta values"
)
lines(y$z, z$en)
Run the code above in your browser using DataLab