Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- set.seed(123)
- ## true values
- mu <- 0.0
- sd <- 1.0
- n <- 10 # sample size
- nbatches <- 1000 # number of samples
- ## Functions that computes the probability that a reference value trueValue
- ## falls into confidence intervals computed by sampling from a normal with
- ## mean center and standard deviation stdev.
- getProbConf <- function(center, stdev, trueValue, alpha, n, nbatches){
- # center : mean of the normal from which is sampled
- # stdev : standard deviation of the normal
- # trueValue: mean of the distribution from which sample was drawn
- # alpha : significance level
- # n : sample size
- # nbatches : number of samples that will be drawn
- q <- qnorm(p=1-alpha/2)
- tmp <- t(replicate(n=nbatches, expr={
- rn <- center+stdev*rnorm(n)
- ci <- mean(rn) + c(-1,1)*q*stdev/sqrt(n)
- findInterval(x=trueValue, vec=ci)
- }))
- return(mean(tmp == 1))
- }
- ## Probability that mu falls into confidence intervals computed from samples
- ## from a normal N(mu,sd). Should be ~0.95
- getProbConf(center=mu, stdev=sd, trueValue=mu,
- alpha=0.05, n=n, nbatches=nbatches)
- ## We draw some observation from the true distribution
- x <- mu+sqrt(sd)*rnorm(n)
- x_bar <- mean(x)
- sd_hat <- sd(x)
- ## Probability that x_bar falls into confidence intervals computed from samples
- ## from a normal N(x_bar,sd_hat). Should be ~0.95
- getProbConf(center=x_bar, stdev=sd_hat, trueValue=x_bar,
- alpha=0.05, n=n, nbatches=nbatches)
- ## Probability that mu falls into confidence intervals computed from samples
- ## from a normal N(x_bar,sd_hat).
- getProbConf(center=x_bar, stdev=sd_hat, trueValue=mu,
- alpha=0.05, n=n, nbatches=nbatches)
- ## Do the same for nrep samples from N(mu, sd). Takes 10 seconds...
- nrep <- 1000
- probs <- replicate(nrep, {
- x <- mu+sqrt(sd)*rnorm(n)
- x_bar <- mean(x)
- sd_hat <- sd(x)
- getProbConf(center=x_bar, stdev=sd_hat, trueValue=mu,
- alpha=0.05, n=n, nbatches=nbatches)
- })
- plot(density(probs))
- ## mean
- mean(probs)
- ## mode
- d <- density(probs)
- d$x[which.max(d$y)]
- ## On average, the true mean (mu) from which a sample with mean x_bar and
- ## standard deviation sd_hat was generated will not lie with a probability of
- ## 0.95 in a 0.95-confidence interval constructed from a sample
- ## from N(x_bar, sd_hat)
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement