r:sampling_distribution
Differences
This shows you the differences between two versions of the page.
Both sides previous revisionPrevious revisionNext revision | Previous revision | ||
r:sampling_distribution [2025/09/10 11:29] – [distribution of sample means] hkimscil | r:sampling_distribution [2025/09/10 20:41] (current) – [PS1. week02] hkimscil | ||
---|---|---|---|
Line 1: | Line 1: | ||
====== PS1. week02 ====== | ====== PS1. week02 ====== | ||
- | |||
< | < | ||
rm(list=ls()) | rm(list=ls()) | ||
Line 32: | Line 31: | ||
var(p1) | var(p1) | ||
- | hist(p1, breaks=100, col=rgb(1, | + | |
+ | hist(p1, breaks=50, col = rgb(1, 1, 1, 0.5), | ||
+ | main = " | ||
+ | abline(v=mean(p1), | ||
+ | hist(p2, add=T, breaks=50, col=rgb(1, | ||
+ | abline(v=mean(p2), | ||
+ | |||
+ | |||
+ | hist(p1, breaks=50, col=rgb(0, | ||
abline(v=mean(p1), | abline(v=mean(p1), | ||
abline(v=mean(p1)-sd(p1), | abline(v=mean(p1)-sd(p1), | ||
Line 90: | Line 97: | ||
pnorm(1.8)-pnorm(-1.8) | pnorm(1.8)-pnorm(-1.8) | ||
- | hist(z.p1, breaks=100, col=rgb(0,0,0,0)) | + | hist(z.p1, breaks=50, col=rgb(1,0,0,0)) |
abline(v=c(m.p1, | abline(v=c(m.p1, | ||
1-(pnorm(1.8)-pnorm(-1.8)) | 1-(pnorm(1.8)-pnorm(-1.8)) | ||
Line 104: | Line 111: | ||
# | # | ||
- | hist(p1, breaks=100, col=rgb(1,1,1,1)) | + | hist(p1, breaks=50, col=rgb(.9,.9,.9,.9)) |
abline(v=mean(p1), | abline(v=mean(p1), | ||
abline(v=mean(p1)-sd(p1), | abline(v=mean(p1)-sd(p1), | ||
Line 116: | Line 123: | ||
c(a, b) | c(a, b) | ||
c(-1, 1) | c(-1, 1) | ||
+ | # note that | ||
+ | .32/2 | ||
+ | pnorm(-1) | ||
+ | qnorm(.32/ | ||
+ | qnorm(pnorm(-1)) | ||
# 95% | # 95% | ||
Line 122: | Line 134: | ||
c(c, d) | c(c, d) | ||
c(-2,2) | c(-2,2) | ||
+ | |||
# 99% | # 99% | ||
e <- qnorm(.01/ | e <- qnorm(.01/ | ||
Line 127: | Line 140: | ||
c(e,f) | c(e,f) | ||
c(-3,3) | c(-3,3) | ||
+ | |||
pnorm(b)-pnorm(a) | pnorm(b)-pnorm(a) | ||
Line 140: | Line 154: | ||
################################ | ################################ | ||
- | hist(p1, breaks=50, col = rgb(1, 0, 0, 0.5), | ||
- | main = " | ||
- | abline(v=mean(p1), | ||
- | hist(p2, add=T, breaks=50, col=rgb(0, | ||
- | abline(v=mean(p2), | ||
- | |||
s.size <- 10 | s.size <- 10 | ||
Line 170: | Line 178: | ||
se.s <- sd(means) | se.s <- sd(means) | ||
- | hist(means, breaks=100, col=rgb(.1, 0, 0, .5)) | + | hist(means, breaks=50, |
- | abline(v=mean(means), | + | xlim = c(mean(means)-5*sd(means), |
+ | col=rgb(1, | ||
+ | abline(v=mean(means), | ||
# now we want to get sd of this distribution | # now we want to get sd of this distribution | ||
lo1 <- mean(means)-se.s | lo1 <- mean(means)-se.s | ||
Line 180: | Line 189: | ||
lo3 <- mean(means)-3*se.s | lo3 <- mean(means)-3*se.s | ||
hi3 <- mean(means)+3*se.s | hi3 <- mean(means)+3*se.s | ||
- | |||
- | hist(means, | ||
- | xlim = c(mean(means)-5*sd(means), | ||
- | col = rgb(1, 0, 0, .5)) | ||
abline(v=mean(means), | abline(v=mean(means), | ||
- | # abline(v=mean(p2), | + | # abline(v=mean(p2), |
abline(v=c(lo1, | abline(v=c(lo1, | ||
- | | + | |
| | ||
Line 198: | Line 203: | ||
# sd of sample means (sd(means)) | # sd of sample means (sd(means)) | ||
- | # is sqrt(var(s1)/ | ||
- | # sd(s1) / sqrt(s.size) | ||
# = se.s | # = se.s | ||
# when iter value goes to | # when iter value goes to | ||
- | # unlimited | + | # infinite |
# mean(means) = mean(p1) | # mean(means) = mean(p1) | ||
# and | # and | ||
# sd(means) = sd(p1) / sqrt(s.size) | # sd(means) = sd(p1) / sqrt(s.size) | ||
- | # that is, sd(means) | + | # that is, se.s = se.z |
# This is called CLT (Central Limit Theorem) | # This is called CLT (Central Limit Theorem) | ||
+ | # see http:// | ||
+ | |||
mean(means) | mean(means) | ||
mean(p1) | mean(p1) | ||
sd(means) | sd(means) | ||
var(p1) | var(p1) | ||
+ | # remember we started talking sample size 10 | ||
sqrt(var(p1)/ | sqrt(var(p1)/ | ||
se.z | se.z | ||
Line 237: | Line 243: | ||
- | hist(means, | + | hist(means, breaks=50, |
xlim = c(mean(means)-5*sd(means), | xlim = c(mean(means)-5*sd(means), | ||
- | col = rgb(1, | + | col = rgb(1, |
- | abline(v=mean(means), | + | abline(v=mean(means), |
# abline(v=mean(p2), | # abline(v=mean(p2), | ||
abline(v=c(lo1, | abline(v=c(lo1, | ||
- | | + | |
| | ||
Line 257: | Line 263: | ||
m.sample.i.got | m.sample.i.got | ||
- | hist(means, | + | hist(means, breaks=30, |
- | xlim = c(mean(means)-10*sd(means), mean(means)+10*sd(means)), | + | xlim = c(mean(means)-7*sd(means), mean(means)+10*sd(means)), |
- | col = rgb(1, | + | col = rgb(1, |
abline(v=mean(means), | abline(v=mean(means), | ||
abline(v=m.sample.i.got, | abline(v=m.sample.i.got, | ||
Line 276: | Line 282: | ||
# (green line) | # (green line) | ||
tmp <- mean(means) - (m.sample.i.got - mean(means)) | tmp <- mean(means) - (m.sample.i.got - mean(means)) | ||
- | abline(v=tmp, | + | abline(v=tmp, |
2 * pnorm(m.sample.i.got, | 2 * pnorm(m.sample.i.got, | ||
m.sample.i.got | m.sample.i.got | ||
### one more time | ### one more time | ||
+ | # this time, with a story | ||
mean(p2) | mean(p2) | ||
sd(p2) | sd(p2) | ||
Line 287: | Line 294: | ||
m.sample.i.got | m.sample.i.got | ||
- | hist(means, | + | tmp <- mean(means) - (m.sample.i.got-mean(means)) |
- | xlim = c(mean(means)-15*sd(means), | + | tmp |
- | col = rgb(1, | + | |
- | abline(v=mean(means), | + | hist(means, breaks=30, |
- | abline(v=m.sample.i.got, | + | xlim = c(tmp-4*sd(means), |
+ | col = rgb(1, | ||
+ | abline(v=mean(means), | ||
+ | abline(v=m.sample.i.got, | ||
# what is the probablity of getting | # what is the probablity of getting | ||
Line 304: | Line 314: | ||
# mean(means) - m.sample.i.got - mean(means) | # mean(means) - m.sample.i.got - mean(means) | ||
# (green line) | # (green line) | ||
- | tmp <- mean(means) - (m.sample.i.got - mean(means)) | + | abline(v=tmp, |
- | abline(v=tmp, | + | |
2 * pnorm(m.sample.i.got, | 2 * pnorm(m.sample.i.got, | ||
- | |||
- | |||
- | |||
- | |||
</ | </ | ||
====== output ====== | ====== output ====== | ||
Line 644: | Line 649: | ||
</ | </ | ||
<WRAP column half> | <WRAP column half> | ||
- | .... | + | qnorm |
* qnorm는 pnorm의 반대값을 구하는 명령어 | * qnorm는 pnorm의 반대값을 구하는 명령어 | ||
* 히스토그램에서 검정 색 부분의 바깥 쪽 부분은 32%이고 왼 쪽의 것은 이것의 반인 16% 이다. | * 히스토그램에서 검정 색 부분의 바깥 쪽 부분은 32%이고 왼 쪽의 것은 이것의 반인 16% 이다. | ||
Line 871: | Line 876: | ||
[1] 3.162278 | [1] 3.162278 | ||
</ | </ | ||
- | * 사실, 우리가 백만 번의 샘플을 취해서 구한 means 집합의 평균과 표준편차 값은 | + | * 사실, 우리가 백만 번의 샘플을 취해서 구한 means 집합의 평균과 표준편차 값은 |
+ | * 이것을 정확하게는 mean of the distribution of sample means, | ||
+ | * standard deviation of the distribution of sample means 라고 부른다 | ||
+ | * 그리고 standard deviation of the distribution of sample means를 흔히 | ||
+ | * standard error라고 부른다. | ||
+ | * script에서 se.s 그리고 se.z 와 같은 se를 쓴 이유도 standard error 값을 구한 것이기 때문이다. | ||
* 만약에 백만 번이 아니라 무한 대로 더 큰 숫자를 사용한다고 하면 | * 만약에 백만 번이 아니라 무한 대로 더 큰 숫자를 사용한다고 하면 | ||
* 위의 se.z 값을 구하는 식의 값을 갖게 된다. 이것을 말로 풀어서 설명하면 | * 위의 se.z 값을 구하는 식의 값을 갖게 된다. 이것을 말로 풀어서 설명하면 | ||
* 샘플평균들의 집합에서 표준편차 값은 | * 샘플평균들의 집합에서 표준편차 값은 | ||
- | * 원래 모집단의 분산값을 샘플사이즈로 나누어준 값에 제곱근을 씌워서 구할 수 있다이다. | + | * 원래 모집단의 분산값을 샘플사이즈로 나누어준 값에 제곱근을 씌워서 구할 수 있다이다 |
+ | |||
+ | \begin{eqnarray*} | ||
+ | \overline{X} \sim \displaystyle \text{N} \left(\mu, \dfrac{\sigma^{2}}{n} \right) | ||
+ | \end{eqnarray*} | ||
+ | |||
+ | \begin{eqnarray} | ||
+ | & & \text{Normal distribution of sample means.} \\ | ||
+ | & & \mu_{\overline{X}} = \mu \\ | ||
+ | & & (\sigma_{\overline{X}})^2 = \frac{\sigma^2}{n} \;\; \text{or} \;\; \sigma_{\overline{X}} = \frac{\sigma}{\sqrt{n}} | ||
+ | \end{eqnarray} | ||
+ | |||
+ | see [[:Central Limit Theorem]] and [[:Sampling Distribution]] | ||
* 즉, 샘플평균을 모은 집합의 분산값은 그 샘플이 추출된 원래 population의 분산값을 샘플크기로 (sample size) 나누어 준 값이다. | * 즉, 샘플평균을 모은 집합의 분산값은 그 샘플이 추출된 원래 population의 분산값을 샘플크기로 (sample size) 나누어 준 값이다. | ||
Line 999: | Line 1021: | ||
</ | </ | ||
+ | ===== Last one . . . Important ===== | ||
<WRAP group> | <WRAP group> | ||
<WRAP column half> | <WRAP column half> | ||
Line 1009: | Line 1031: | ||
[1] 10 | [1] 10 | ||
> sample(p2, s.size) | > sample(p2, s.size) | ||
- | | + | |
+ | [9] 116.41683 105.91887 | ||
> m.sample.i.got <- mean(sample(p2, | > m.sample.i.got <- mean(sample(p2, | ||
> m.sample.i.got | > m.sample.i.got | ||
- | [1] 120.2492 | + | [1] 114.155 |
- | > | + | > temp <- mean(means) - (m.sample.i.got-mean(means)) |
- | > hist(means, | + | > temp |
- | + xlim = c(mean(means)-15*sd(means), | + | [1] 85.84835 |
+ | > hist(means, breaks=30, | ||
+ | + xlim = c(temp-4*sd(means), | ||
+ col = rgb(1, 0, 0, .5)) | + col = rgb(1, 0, 0, .5)) | ||
> abline(v=mean(means), | > abline(v=mean(means), | ||
- | > abline(v=m.sample.i.got, | + | > abline(v=m.sample.i.got, |
- | > | + | > abline(v=temp, |
> # what is the probablity of getting | > # what is the probablity of getting | ||
> # greater than | > # greater than | ||
> # m.sample.i.got? | > # m.sample.i.got? | ||
> m.sample.i.got | > m.sample.i.got | ||
- | [1] 120.2492 | + | [1] 114.155 |
> pnorm(m.sample.i.got, | > pnorm(m.sample.i.got, | ||
- | [1] 7.560352e-11 | + | [1] 3.820806e-06 |
- | > | + | |
> # then, what is the probabilty of getting | > # then, what is the probabilty of getting | ||
> # greater than m.sample.i.got and | > # greater than m.sample.i.got and | ||
Line 1036: | Line 1060: | ||
> abline(v=tmp, | > abline(v=tmp, | ||
> 2 * pnorm(m.sample.i.got, | > 2 * pnorm(m.sample.i.got, | ||
- | [1] 1.51207e-10 | + | [1] 7.641611e-06 |
> | > | ||
</ | </ | ||
Line 1042: | Line 1066: | ||
<WRAP column half> | <WRAP column half> | ||
.... | .... | ||
+ | 다른 모집단에서 온 샘플 (sample from other population, p2) | ||
+ | 이전 처럼 10명을 샘플로 뽑는데, p2에서 뽑는다. 따라서 이 샘플의 평균은 114.155 이다. 그런데 연구자는 이 샘플이 나온 모집단은 전혀 모르는 상태이다. 모집단의 평균도 모르고 분산값도 모른다. 그냥 10명을 샘플로 뽑았고, 이들이 p1에서 나왔을 확률을 알아보려고 한다. 이를 알아보기 위해서 스크립트를 돌려 보니, blue와 green라인 밖의 범위가 나올 확률은 7.641611e-06로 0에 가깝다. 따라서, 이 샘플은 10명은 p1에 속한 원소가 아닌 다른 성격을 가진 원소이다. | ||
+ | |||
+ | |||
+ | 위에서 p2의 parameter에 대해서 잘 모른다는 점에 주목하라. 그리고 아래 시나리오를 상상하라. | ||
+ | 어느 한 모집단의 IQ 평균이 100 이고 표준편차가 10 임을 알고 있다. 확률과통계 교수는 머리가 좋아지는 약을 개발하여 이를 팔아보려고 하고 있다. 이를 위해서 확통교수는 10을 뽑아서 머리가 좋아지는 약을 복용하도록 한 후에 IQ를 측정하였다. 그런데, 그 IQ 평균이 114.155 이다. 이 점수를 가지고 약의 효과가 있는지 검증을 해보력고 한다. | ||
</ | </ | ||
</ | </ | ||
+ | {{: | ||
r/sampling_distribution.1757471378.txt.gz · Last modified: 2025/09/10 11:29 by hkimscil