-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathweek_05_assign_05.Rmd
147 lines (106 loc) · 3.67 KB
/
week_05_assign_05.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
---
title: "week_05_assign_05"
author: "[email protected]"
date: "`r Sys.Date()`"
output: html_document
---
let us first compute the probability corresponding to the below quantile using the *pnorm function* and the given normal distribution parameters
```{r}
# P(X < 12)
quan_12 <- pnorm(q = 12, mean = 10, sd = sqrt(4))
quan_12
# P(X > 10)
quan_10 <- pnorm(q = 10 , mean =10 , sd = sqrt(4))
quan_10
#P( 10 < X < 12) = F(12) - F(10)
quan_12-quan_10
```
let us compute the 90th percentile of an expo distribution with rate parameter 1/2 , 1/4 and 2
```{r}
qexp(p = 0.9, rate = 1/2)
qexp(p = 0.9, rate = 1/4)
qexp(p = 0.9, rate = 2)
```
```{r}
# P(X=2)
dpois(x = 2,lambda = 4)
# P(X <= 2)
# Careful with the definition of CDF: For <= 2, we use ppois() or we add P(X=0)+P(X=1)+P(X=2)
ppois(q = 2,lambda = 4, lower.tail = TRUE)
#dpois(x = 0, lambda = 4) + dpois(x = 1, lambda = 4) + dpois(x = 2, lambda = 4) alternatively
sum(dpois(x = 0:2,lambda = 4))
```
let us generate some random numbers from normal with parameters
1.5 and 2 as mean and standard deviation respectively
```{r}
norm.sample <- rnorm(n = 200, mean = 1.5, sd = 2)
head(norm.sample)
max(norm.sample);min(norm.sample);range(norm.sample);norm.sample[76:79]
```
```{r}
plot(density(norm.sample), main = "Kernel Density")
```
```{r}
plot(density(norm.sample), main = "Kernel Density")
lines(norm.sample, dnorm(norm.sample, mean =1.5, sd = 2 ), col="red")
```
let us do some ordering or sorting the points
```{r}
plot(density(norm.sample), main = "Kernel Density")
sorted.norm.sample = sort(norm.sample)
lines(sorted.norm.sample, dnorm(sorted.norm.sample,mean = 1.5, sd = 2),
col="red")
```
note that above the black line represents the simulated data points and the red is the one which represents actually the theoretical normal distribution
```{r}
mean(sorted.norm.sample);sd(sorted.norm.sample)
mean(norm.sample);sd(norm.sample)
```
As an example, we shall study the behavior of the sampling distribution of the median of a sample, median $\left(X_1, \ldots, X_n\right)$, when the observations are drawn from an exponential distribution, with the below CDF
$$
F(x)=1-\exp (-\lambda x)
$$
To look at the sampling distribution of the median when $\lambda=1 / 2$ and sample size is 100, we simulate 500 samples: we proceed via standard for loop
```{r}
exp.median <- numeric(500)
for (i in 1:500){
exp.median[i]<- median(rexp(100, rate = 1/2))
}
```
Explore graphically the sampling distribution of the above simulated median (with histograms, boxplots etc).
```{r}
par(mfrow=c(1,3))
hist(exp.median, probability = TRUE, col = "darkseagreen")
#lines(density(Exp.median),col="navyblue")
boxplot(exp.median, probability = TRUE, col = "salmon")
abline(h=1.39, col = "navyblue") #theoretical median we computed in point 1.
library(vioplot)
vioplot(exp.median, col = "lightblue")
```
let us check normality Q-Q plot
```{r}
par(mfrow=c(1,1))
qqnorm(exp.median)
qqline(exp.median)
# certainly looks like normal
```
Redoing the same steps as before for the 90-th quantile of the distribution rather than the 50-th percentile
```{r}
exp_90q <- numeric(500)
for (i in 1:500){
exp_90q [i]<- quantile(rexp(100, rate = 1/2),probs = 0.9)
} # where number of experiments = 500 and sample size per experiment is 100
```
```{r}
par(mfrow=c(1,3))
hist(exp_90q, probability = TRUE, col = "darkseagreen")
boxplot(exp_90q, probability = TRUE, col = "salmon")
abline(h=1.39, col = "navyblue")
library(vioplot)
vioplot(exp_90q, col = "lightblue")
```
```{r}
par(mfrow=c(1,1))
qqnorm(exp_90q)
qqline(exp_90q)
```