-
Notifications
You must be signed in to change notification settings - Fork 0
/
Monte_Carlo_Confidence_Interval_Normal_Distribution.rmd
265 lines (187 loc) · 8.18 KB
/
Monte_Carlo_Confidence_Interval_Normal_Distribution.rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
---
title: "Using Monte Carlo simulation to estimate the confidence interval and estimate the length of the confidence interval in a normal population in the cases of equal and unequal tails, as well as the standard deviation being known and the standard deviation being unknown"
author: "habib Ezzatabadi"
output: github_document
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```
## Define a function for simulate CI and get length of CI in Normal population with equal tails and unequal tails
```{r}
lfun_ci <- function(
mu = 1, ## Mean of population
sig = 1, ## standard devitaion of population
state_sig = NA,
## if state_sig = NULL, function calculate CI with
#unknown statndard deviation
nsim = 1e+4, ## size of simulation, default = 10000
nsamp = 20, ## size of sample, default = 20
alpha = 0.05, ## alpha value for get 100(1-alpha) % oc CI,
## alpha default = 0.05
left_tail = 1/2,
seed = 12334 # for get same results
) {
set.seed(seed)
Mat_sim <- matrix(rnorm(nsim * nsamp, mean = mu, sd = sig),
nsim, nsamp) ## simulate matrix of samples
## calculate standard deviation of simulated samples
ss <- apply(Mat_sim, 1, sd)
## calculate mean of simulated samples
Xbar <- rowMeans(Mat_sim)
a1 <- alpha * left_tail ## calculate alpha for lower tail
a2 <- alpha - a1 ## calculate alpha for upper tail
## get z quantile for left bond
z1 <- qnorm(a1, mean = 0, sd = 1, lower.tail = FALSE)
## get z quantile for right bond
z2 <- qnorm(a2, , mean = 0, sd = 1, lower.tail = FALSE)
## get t quantile for left bond
t1 <- qt(a1, df = nsamp-1, lower.tail = FALSE)
## get t quantile for right bond
t2 <- qt(a2, df = nsamp-1, lower.tail = FALSE)
if (! is.null(state_sig)) { ## check standard
## deviatioan known or unknown
## get lower bond of CI
lower_bond = Xbar - z1 * sig/sqrt(nsamp)
## get upper bond of CI
upper_bond = Xbar + z2 * sig/sqrt(nsamp)
} else { ## if standard deviation is unknown use these bonds
## get lower bond
lower_bond = Xbar - t1 * ss / sqrt(nsamp)
## get upper bond
upper_bond = Xbar + t2 * ss / sqrt(nsamp)
}
L <- upper_bond - lower_bond ## calculate length of bonds
## Calculate the confidence percentage of the
## simulated confidence interval
## condition 1: mu > lower_bond
temp1 <- 1 * (mu > lower_bond)
## condition 2: mu < upper_bond
temp2 <- 1 * (mu < upper_bond)
## get intersection between condition 1 and condition 2
temp3 <- temp1 * temp2
ci_sim <- mean(temp3)
## calculate output of function
Result <- list(
Tails = sprintf("Tails : %s", ifelse(left_tail == 1/2,
"Equal", "Unequal")),
alpha = sprintf("alpha = %.2f", alpha),
sample_size = sprintf("size of sample = %d", nsamp),
simulation_size =
sprintf("size of Simulation = %d", nsim),
standard_deviation =
sprintf("Standard Deviation = %s",
ifelse(is.null(state_sig), "Unknown",
as.character(sig))),
Mean = sprintf("Mean = %.2f", mu),
Result = c("Average of CI simulation" = ci_sim,
"Average of Length of CI" = mean(L)))
return (Result)
}
```
***
***
## Example I:
#### $\mu = 2~~\sigma = 3$, sample size = 15, sample of simulation = 10000,
#### Standard Deviation = Known and Tails are Equal
```{r}
### EX I
## Mean = 2, Standard Deviation = 3, sample size = 15
lfun_ci(nsamp = 15, mu = 2, sig = 3)
```
***
## Example II:
#### $\mu = 10, ~~\sigma = 2.5$, sample size = 15,
#### size of simulation = 1000, Standard Deviaion = Known and Tails are Equal
```{r}
lfun_ci(nsamp = 10, mu = 10, sig = 2.5, nsim = 1000)
```
***
## Example III:
#### $\mu = 2, ~~\sigma = 3$, sample size = 15,
#### size of simulation = 10000, Standard Deviaion = Known and tails are Unequal
```{r}
## Mean = 2, Standard Deviation = 3, sample size = 15,
## alpha of left tail = alpha / 3
lfun_ci(nsamp = 15, mu = 2, sig = 3, left_tail = 1 / 3)
```
***
## Example IV:
#### $\mu = 10, ~~\sigma = 2.5$, sample size = 10,
#### size of simulation = 1000, Standard Deviaion = Known and tails are Unequal
```{r}
## Mean = 10, Standard Deviation = 2.5, sample size = 15,
## size of simulation = 1000, alpha of left tail = alpha / 4
lfun_ci(nsamp = 10, mu = 10, sig = 2.5, nsim = 1000,
left_tail = 1/4)
```
***
## Example V:
#### $\mu = 2, ~~\sigma = 3$, sample size = 15,
#### size of simulation = 10000, Standard Deviaion = UnKnown and tails are Equal
```{r}
## Mean = 2, Standard Deviation = 3, sample size = 15,
## Standard Deviation is Unknown
lfun_ci(nsamp = 15, mu = 2, sig = 3, state_sig = NULL)
```
***
## Example VI:
#### $\mu = 10, ~~\sigma = 2.5$, sample size = 10,
#### size of simulation = 1000, Standard Deviaion = UnKnown and tails are Equal
```{r}
## Mean = 10, Standard Deviation = 2.5, sample size = 15,
## size of simulation = 1000, Standard Deviation is Unknown
lfun_ci(nsamp = 10, mu = 10, sig = 2.5, nsim = 1000, state_sig = NULL)
```
***
## Example VII:
#### $\mu = 2, ~~\sigma = 3$, sample size = 15,
#### size of simulation = 10000, Standard Deviaion = UnKnown and tails are Unequal
```{r}
## Mean = 2, Standard Deviation = 3, sample size = 15,
## alpha of left tail = alpha / 3, Standard Deviaion is Unknown
lfun_ci(nsamp = 15, mu = 2, sig = 3, left_tail = 1 / 3, state_sig = NULL)
```
***
## Example VIII:
#### $\mu = 10, ~~\sigma = 2.5$, sample size = 10,
#### size of simulation = 1000, Standard Deviaion = UnKnown and tails are Unequal
```{r}
## Mean = 10, Standard Deviation = 2.5, sample size = 15,
## size of simulation = 1000, alpha of left tail = alpha / 4,
## Standard Deviation is Unknown
lfun_ci(nsamp = 10, mu = 10, sig = 2.5, nsim = 1000,
left_tail = 1/4, state_sig = NULL)
```
***
***
## Compare results
```{r}
## compare results
temp1 <- lfun_ci(nsamp = 15, mu = 2, sig = 3)$Result
temp2 <- lfun_ci(nsamp = 15, mu = 2, sig = 3, left_tail = 1 / 3)$Result
temp3 <- lfun_ci(nsamp = 15, mu = 2, sig = 3, state_sig = NULL)$Result
temp4 <- lfun_ci(nsamp = 15, mu = 2, sig = 3, left_tail = 1 / 3,
state_sig = NULL)$Result
Compare_Result_1 <- rbind(temp1, temp2, temp3, temp4)
row_names <- c("Known sd with equal tails",
"Known sd with unequal tails",
"UnKnown sd with equal tails",
"UnKnown sd with Unequal tails")
res1 <- as.data.frame(Compare_Result_1)
rownames(res1) <- row_names
knitr :: kable(res1, align = "c", caption = "Compare of Results for
sample size = 15, Mean = 2, Standard Deviation = 3")
## For Ex II
temp21 <- lfun_ci(nsamp = 10, mu = 10, sig = 2.5, nsim = 1000)$Result
temp22 <- lfun_ci(nsamp = 10, mu = 10, sig = 2.5, nsim = 1000,
left_tail = 1/4)$Result
temp23 <- lfun_ci(nsamp = 10, mu = 10, sig = 2.5, nsim = 1000,
state_sig = NULL)$Result
temp24 <- lfun_ci(nsamp = 10, mu = 10, sig = 2.5, nsim = 1000,
left_tail = 1/4, state_sig = NULL)$Result
Compare_Result_2 <- rbind(temp21, temp22, temp23, temp24)
res2 <- as.data.frame(Compare_Result_2)
rownames(res2) <- row_names
knitr :: kable(res2, align = "c", caption = "Compare of Results for
sample size = 10, Mean = 10, Standard Deviation = 2.5")
```