• API
• FAQ
• Tools
• Archive
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.
Top