SHARE
TWEET

Untitled

a guest Nov 28th, 2014 309 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. set.seed(123)
  2. ## true values
  3. mu <- 0.0
  4. sd <- 1.0
  5.  
  6. n <- 10 # sample size
  7. nbatches <- 1000 # number of samples
  8.  
  9.  
  10. ## Functions that computes the probability that a reference value trueValue
  11. ## falls into confidence intervals computed by sampling from a normal with
  12. ## mean center and standard deviation stdev.
  13. getProbConf <- function(center, stdev, trueValue, alpha, n, nbatches){
  14.   # center   : mean of the normal from which is sampled
  15.   # stdev    : standard deviation of the normal
  16.   # trueValue: mean of the distribution from which sample was drawn
  17.   # alpha    : significance level
  18.   # n        : sample size
  19.   # nbatches : number of samples that will be drawn  
  20.   q <- qnorm(p=1-alpha/2)
  21.   tmp <- t(replicate(n=nbatches, expr={
  22.     rn <- center+stdev*rnorm(n)
  23.     ci <- mean(rn) + c(-1,1)*q*stdev/sqrt(n)
  24.     findInterval(x=trueValue, vec=ci)
  25.   }))
  26.   return(mean(tmp == 1))
  27. }
  28.  
  29.  
  30. ## Probability that mu falls into confidence intervals computed from samples
  31. ## from a normal N(mu,sd). Should be ~0.95
  32. getProbConf(center=mu, stdev=sd, trueValue=mu,
  33.             alpha=0.05, n=n, nbatches=nbatches)
  34.  
  35.  
  36. ## We draw some observation from the true distribution
  37. x <- mu+sqrt(sd)*rnorm(n)
  38. x_bar <- mean(x)
  39. sd_hat <- sd(x)
  40.  
  41.  
  42. ## Probability that x_bar falls into confidence intervals computed from samples
  43. ## from a normal N(x_bar,sd_hat). Should be ~0.95
  44. getProbConf(center=x_bar, stdev=sd_hat, trueValue=x_bar,
  45.             alpha=0.05, n=n, nbatches=nbatches)
  46.  
  47.  
  48.  
  49.  
  50. ## Probability that mu falls into confidence intervals computed from samples
  51. ## from a normal N(x_bar,sd_hat).
  52. getProbConf(center=x_bar, stdev=sd_hat, trueValue=mu,
  53.             alpha=0.05, n=n, nbatches=nbatches)
  54.  
  55.  
  56. ## Do the same for nrep samples from N(mu, sd). Takes 10 seconds...
  57. nrep <- 1000
  58. probs <- replicate(nrep, {
  59.   x <- mu+sqrt(sd)*rnorm(n)
  60.   x_bar <- mean(x)
  61.   sd_hat <- sd(x)
  62.   getProbConf(center=x_bar, stdev=sd_hat, trueValue=mu,
  63.               alpha=0.05, n=n, nbatches=nbatches)
  64. })
  65.  
  66. plot(density(probs))
  67. ## mean
  68. mean(probs)
  69. ## mode
  70. d <- density(probs)
  71. d$x[which.max(d$y)]
  72.  
  73. ## On average, the true mean (mu) from which a sample with mean x_bar and
  74. ## standard deviation sd_hat was generated will not lie with a probability of
  75. ## 0.95 in a 0.95-confidence interval constructed from a sample
  76. ## from N(x_bar, sd_hat)
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
Top