COMMunication
RESearch.NET

User Tools

Site Tools

You are here: index » r » sampling_distribution
Trace:
r:sampling_distribution

Differences

This shows you the differences between two versions of the page.

Link to this comparison view

Both sides previous revisionPrevious revision
Next revision		Previous revision
r:sampling_distribution [2026/03/21 10:28] – [R script output] hkimscilr:sampling_distribution [2026/03/24 10:17] (current) – hkimscil
Line 62: Line 62:
 </WRAP> </WRAP>
 </WRAP> </WRAP>
-pnorm +===== pnorm ===== 
 + 
 <WRAP group> <WRAP group>
 <WRAP column half> <WRAP column half>
Line 182: Line 183:
 </WRAP> </WRAP>
 </WRAP> </WRAP>
-z score, 표준점수 +===== z score, 표준점수 ===== 
 <WRAP group> <WRAP group>
 <WRAP column half> <WRAP column half>
Line 260: Line 262:
  
 ---- ----
-qnorm + 
 +===== qnorm =====
 <WRAP group> <WRAP group>
 <WRAP column half> <WRAP column half>
Line 369: Line 372:
 </WRAP> </WRAP>
 </WRAP> </WRAP>
- 
 ---- ----
-distribution of sample means + 
 +<code> 
 +pnorm(110, 100, 10, lower.tail = F)  
 +pnorm((110-100)/10, lower.tail = F) 
 +pnorm(1, lower.tail = F) 
 +1-pnorm(1) 
 +pnorm(-1) 
 + 
 +1-(pnorm(-1)*2) 
 +</code> 
 +===== distribution of sample means ===== 
 + 
 아래는 모두 같은 의미이다.  아래는 모두 같은 의미이다. 
   * distribution of sample means   * distribution of sample means
Line 465: Line 478:
 <code> <code>
 > mean(means) > mean(means)
-[1] 99.99912+[1] 100.0025
 > sd(means) > sd(means)
-[1] 3.164792+[1] 3.167479
 > var(means) > var(means)
-[1] 10.01591+[1] 10.03292
 </code> </code>
   * 위와 같다. 이 값은 샘플평균들의 평균이 무엇인가와 (mean(means))   * 위와 같다. 이 값은 샘플평균들의 평균이 무엇인가와 (mean(means))
Line 556: Line 569:
 > # that is, se.s = se.z > # that is, se.s = se.z
 > # This is called CLT (Central Limit Theorem) > # This is called CLT (Central Limit Theorem)
-> # see http://commres.net/wiki/cetral_limit_theorem+> # see http://commres.net/central_limit_theorem
 >  > 
 > mean(means) > mean(means)
Line 604: Line 617:
 [[:mean and variance of the sample mean]] [[:mean and variance of the sample mean]]
 .... ....
-  * 그런데 이 값은 (se.s = 3.161886) +  * 샘플들을 직접 모아서 구한 샘플평균들의 표준편차 값은 (standard deviation of sample means, se.s = 3.161886) 
-  * se.z 를 구하는 방법과 거의 같은 값을 갖는다 3.162278+  * se.z 를 구해서 얻은 값과 같다 (3.162278).
 <code> <code>
 > se.z <- sqrt(var(p1)/s.size) > se.z <- sqrt(var(p1)/s.size)
Line 663: Line 676:
 <WRAP column half> <WRAP column half>
 <code> <code>
-> hist(means,  + 
-+      xlim = c(mean(means)-5*sd(means), mean(means)+10*sd(means)),  +> means.10 <- rnorm2(iter, mean(p1), c(se.z)) 
-+      col = rgb(1, 0, 0, .5)) +> hist(means.10, breaks=50, 
-> abline(v=mean(means), col="black", lwd=2)++      xlim = c(mean(p1)-5*se.z, mean(p1)+10*se.z),  
 ++      col = rgb(1, 1, 1, .5)) 
 +> abline(v=mean(p1), col="black", lwd=3)
 > # abline(v=mean(p2), colo='darkgreen', lwd=3) > # abline(v=mean(p2), colo='darkgreen', lwd=3)
-> abline(v=c(lo1, hi1, lo2, hi2, lo3, hi3),  +> abline(v=c(loz1, hiz1, loz2, hiz2, loz3, hiz3),  
-+        col=c("green","green", "blue", "blue", "orange", "orange"), ++        col=c("darkgreen","darkgreen", "blue", "blue", "orange", "orange"), 
 +        lwd=2) +        lwd=2)
 >  > 
Line 686: Line 701:
 [1]  91 109 [1]  91 109
 >  > 
 +>
 +>
 </code> </code>
 </WRAP> </WRAP>
Line 693: Line 710:
 </WRAP> </WRAP>
 </WRAP> </WRAP>
 +
 +===== Hypothesis test =====
  
 <WRAP group> <WRAP group>
 <WRAP column half> <WRAP column half>
 <code> <code>
->  
 > m.sample.i.got <- mean(means)+ 1.5*sd(means) > m.sample.i.got <- mean(means)+ 1.5*sd(means)
 > m.sample.i.got > m.sample.i.got
-[1] 104.7537+[1] 104.725
 >  > 
-> hist(means, breaks=30,  +> hist(means.10, breaks=30,  
-+      xlim = c(mean(means)-7*sd(means), mean(means)+10*sd(means)), ++      xlim = c(mean(p1)-7*se.z, mean(p1)+10*se.z), 
 +      col = rgb(1, 1, 1, .5)) +      col = rgb(1, 1, 1, .5))
-> abline(v=mean(means), col="black", lwd=3)+> abline(v=mean(p1), col="black", lwd=3)
 > abline(v=m.sample.i.got, col='darkgreen', lwd=3) > abline(v=m.sample.i.got, col='darkgreen', lwd=3)
 >  > 
Line 712: Line 730:
 > # m.sample.i.got? > # m.sample.i.got?
 > m.sample.i.got > m.sample.i.got
-[1] 104.7537 +[1] 104.725
-> pnorm(m.sample.i.got, mean(means), sd(means), lower.tail = F) +
-[1] 0.0668072+
 > pnorm(m.sample.i.got, mean(p1), se.z, lower.tail = F) > pnorm(m.sample.i.got, mean(p1), se.z, lower.tail = F)
-[1] 0.06638638+[1] 0.06756624
 >  > 
 > # then, what is the probabilty of getting  > # then, what is the probabilty of getting 
Line 723: Line 739:
 > # 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))+> tmp <- mean(p1) - (m.sample.i.got - mean(p1)) 
 +> tmp 
 +[1] 95.27504
 > abline(v=tmp, col='red', lwd=3) > abline(v=tmp, col='red', lwd=3)
-> 2 * pnorm(m.sample.i.got, mean(p1), se.z, lower.tail = F) 
-[1] 0.1327728 
 > m.sample.i.got > m.sample.i.got
-[1] 104.7537 +[1] 104.725 
->+> # 붉은 선의 왼쪽 부분을 구해서 2를 곱한 값의 범위 
 +> prob.of.m.sample.i.got <- 2 * pnorm(m.sample.i.got, mean(p1), se.z, lower.tail = F) 
 +> prob.of.m.sample.i.got 
 +[1] 0.1351325
 >  > 
-> 2 * pnorm(m.sample.i.got, mean(p1), se.z, lower.tail = F) 
-[1] 0.1327728 
 > (m.sample.i.got - mean(p1))/se.z > (m.sample.i.got - mean(p1))/se.z
          [,1]          [,1]
-[1,] 1.503257+[1,] 1.494165
 > z.score <- (m.sample.i.got - mean(p1))/se.z > z.score <- (m.sample.i.got - mean(p1))/se.z
 > pnorm(z.score, lower.tail = FALSE) > pnorm(z.score, lower.tail = FALSE)
            [,1]            [,1]
-[1,] 0.06638638+[1,] 0.06756624
 > 2 * pnorm(z.score, lower.tail = FALSE) > 2 * pnorm(z.score, lower.tail = FALSE)
           [,1]           [,1]
-[1,] 0.1327728+[1,] 0.1351325 
 +>  
 +> # 즉 n = ssize 크기의 샘플을 구했을 때  
 +> # 104.7489 점 오른쪽 그리고, 
 +> # 이에 대응하는 95.27504 지점 왼쪽의  
 +> # 평균이 나오는 probability는  
 +> # 0.1351325 
 >  > 
- 
 </code> </code>
 </WRAP> </WRAP>
Line 750: Line 772:
 .... ....
   * 만약에 내가 한 샘플을 취해서 평균값을 살펴보니    * 만약에 내가 한 샘플을 취해서 평균값을 살펴보니 
-  * m.sample.i.got 값이었다고 하자 (104.7383).+  * m.sample.i.got 값이었다고 하자 (104.725).
   * 이 값보다 큰 값이거나 아니면    * 이 값보다 큰 값이거나 아니면 
   * 이 값에 해당하는 평균 반대편 값보다 작은 값이 값이    * 이 값에 해당하는 평균 반대편 값보다 작은 값이 값이 
   * 나올 확률은 무엇인가?    * 나올 확률은 무엇인가? 
   * 즉, 녹색선과 연두색 선 바깥 쪽 부분의 probability 값은?   * 즉, 녹색선과 연두색 선 바깥 쪽 부분의 probability 값은?
-  * 아래처럼 구해서 13.4% 정도가 된다+  * 아래처럼 구해서 13.5% 정도가 된다
   * 즉, 모집단 p1에서    * 즉, 모집단 p1에서 
   * 무작위로 샘플링을 하여 (see [[:sampling]] and probability sampling)   * 무작위로 샘플링을 하여 (see [[:sampling]] and probability sampling)
   * s.size의 (10) 샘플을 취했는데 그 샘플의 평균이    * s.size의 (10) 샘플을 취했는데 그 샘플의 평균이 
-  * 104.7383 점보다 크거나 혹은  +  * 104.725 점보다 크거나 혹은  
-  * 95.26505 점보다 작을 확률은 13.4% 이다.  +  * 95.27504 점보다 작을 확률은 13.5% 이다. 
-<code>+
  
 +<code>
 > 2 * pnorm(m.sample.i.got, mean(p1), sd(means), lower.tail = F) > 2 * pnorm(m.sample.i.got, mean(p1), sd(means), lower.tail = F)
 [1] 0.1339882 [1] 0.1339882
 </code> </code>
 +
   * 그런데 위의 pnorm 은 표준점수화 해서 생각할 수 있다.    * 그런데 위의 pnorm 은 표준점수화 해서 생각할 수 있다. 
 +
 <code> <code>
-> 2 * pnorm(m.sample.i.got, mean(p1), sd(means), lower.tail = F) +> (m.sample.i.got - mean(p1))/se.z 
-[1] 0.13371 +         [,1] 
-> (m.sample.i.got - mean(p1))/sd(means) +[1,] 1.494165 
-[1] 1.499631 +> z.score <- (m.sample.i.got - mean(p1))/se.z
-> z.score <- (m.sample.i.got - mean(p1))/sd(means)+
 > pnorm(z.score, lower.tail = FALSE) > pnorm(z.score, lower.tail = FALSE)
-[1] 0.06685502+           [,1] 
 +[1,] 0.06756624
 > 2 * pnorm(z.score, lower.tail = FALSE) > 2 * pnorm(z.score, lower.tail = FALSE)
-[1] 0.13371 +          [,1] 
-> +[1,] 0.1351325
 </code> </code>
   * 위처럼 z score를 구해서 pnorm으로 probability를 보는 것을 z-test 라고 한다.    * 위처럼 z score를 구해서 pnorm으로 probability를 보는 것을 z-test 라고 한다. 
Line 783: Line 807:
 </WRAP> </WRAP>
 ---- ----
-Last one . . . Important +
 <WRAP group> <WRAP group>
 <WRAP column half> <WRAP column half>
 <code> <code>
->  
 > ### one more time  > ### one more time 
 > # this time, with a story > # this time, with a story
 +> # 한편 우리는 특정한 평균과 표준편차를 갖는
 +> # p2를 알고 있다 (110, 10)
 > mean(p2) > mean(p2)
 [1] 110 [1] 110
Line 797: Line 822:
 > sample.i.got.p2 <- sample(p2, s.size) > sample.i.got.p2 <- sample(p2, s.size)
 > sample.i.got.p2 > sample.i.got.p2
- [1] 109.64312 117.28119 106.40247 113.29156 124.23996 117.81039  91.88752 + [1] 106.63828 107.09795 108.08661 119.09268 103.38813 118.30613 111.49542 
- [8] 102.20386  95.59136 101.16378+ [8]  94.25277 115.38045  95.49686
 >  > 
 > m.sample.i.got.p2 <- mean(sample.i.got.p2) > m.sample.i.got.p2 <- mean(sample.i.got.p2)
 > m.sample.i.got.p2 > m.sample.i.got.p2
-[1] 107.9515+[1] 107.9235
 >  > 
-> hist(means, breaks=30,  +> # 위의 샘플 평균이 위에서 얘기했던 샘플평균들의 집합에서  
-+      xlim = c(tmp-10*sd(means), m.sample.i.got+10*sd(means)), +> # 나올 확률은 무엇일까? 
 +> hist(means.10, breaks=30,  
 ++      xlim = c(tmp-10*se.z, m.sample.i.got+10*se.z), 
 +      col = rgb(1, 1, 1, .5)) +      col = rgb(1, 1, 1, .5))
-> abline(v=mean(means), col="black", lwd=3)+> abline(v=mean(p1), col="black", lwd=3)
 > abline(v=m.sample.i.got.p2, col='blue', lwd=3) > abline(v=m.sample.i.got.p2, col='blue', lwd=3)
 >  > 
Line 814: Line 841:
 > # m.sample.i.got? > # m.sample.i.got?
 > m.sample.i.got.p2 > m.sample.i.got.p2
-[1] 107.9515 +[1] 107.9235 
-> pnorm(m.sample.i.got.p2, mean(p1), se.z, lower.tail = F) +> pnorm(m.sample.i.got.p2, mean(p1), se.z, lower.tail = F)  
-[1] 0.005960209+[1] 0.006111512 
 +> # 혹은 
 +> z.cal <- (m.sample.i.got.p2-mean(p1))/se.z 
 +> pnorm(z.cal, lower.tail = F) 
 +            [,1] 
 +[1,] 0.006111512
 >  > 
-> # then, what is the probabilty of getting +> # then, what is the probability of getting 
 > # greater than m.sample.i.got and > # greater than m.sample.i.got and
 > # less than corresponding value, which is > # less than corresponding value, which is
Line 825: Line 857:
 > abline(v=mean(p1)-(m.sample.i.got.p2-mean(p1)), col='red', lwd=3) > abline(v=mean(p1)-(m.sample.i.got.p2-mean(p1)), col='red', lwd=3)
 > 2 * pnorm(m.sample.i.got.p2, mean(p1), se.z, lower.tail = F) > 2 * pnorm(m.sample.i.got.p2, mean(p1), se.z, lower.tail = F)
-[1] 0.01192042 +[1] 0.01222302
->  +
-> z.cal <- (m.sample.i.got.p2 - mean(p1)) / se.z +
-> se.z +
-         [,1] +
-[1,] 3.162278 +
-> z.cal +
-         [,1] +
-[1,] 2.514492 +
->  +
-> pnorm(z.cal, lower.tail = F) +
-            [,1] +
-[1,] 0.005960209+
 > 2 * pnorm(z.cal, lower.tail = F) > 2 * pnorm(z.cal, lower.tail = F)
            [,1]            [,1]
-[1,] 0.01192042+[1,] 0.01222302 
 +>  
 +> # 위는 샘플평균과 모집단 평균값 간의 차이 diff를 
 +> # standard error값으로 즉,  
 +> # sd of sample means (n=s.size) 값으로 나눠 준  
 +> # 표준점수 값을 구하고  
 +> # 이 점수에 대한 probability값을 구해 본것
 >  > 
 > diff <- m.sample.i.got.p2 - mean(p1) > diff <- m.sample.i.got.p2 - mean(p1)
Line 847: Line 873:
 > z.cal > z.cal
          [,1]          [,1]
-[1,] 2.514492 +[1,] 2.505639
-> pnorm(z.cal, lower.tail=F)*2 +
-           [,1] +
-[1,] 0.01192042+
 > prob.z <- pnorm(z.cal, lower.tail=F)*2 > prob.z <- pnorm(z.cal, lower.tail=F)*2
-> +> prob.z 
 +           [,1] 
 +[1,] 0.01222302 
 +> # 위에 대한 해석 . . . . .  
 +> # 
 +> # 
 +> ###################################################
 </code> </code>
 </WRAP> </WRAP>
Line 858: Line 887:
 .... ....
 다른 모집단에서 온 샘플 (sample from other population, p2) 다른 모집단에서 온 샘플 (sample from other population, p2)
-이전 처럼 10명을 샘플로 뽑는데, p2에서 뽑는다. 따라서 이 샘플의 평균은 114.155 이다. 그런데 연구자는 이 샘플이 나온 모집단은 전혀 모르는 상태이다. 모집단의 평균도 모르고 분산값도 모른다. 그냥 10명을 샘플로 뽑았고, 이들이 p1에서 나왔을 확률을 알아보려고 한다. 이를 알아보기 위해서 스크립트를 돌려 보니, blue와 green라인 밖의 범위가 나올 확률은 7.641611e-06로 0에 가깝다. 따라서, 이 샘플은 10명은 p1에 속한 원소가 아닌 다른 성격을 가진 원소이다. +이전 처럼 10명을 샘플로 뽑는데, p2에서 뽑는다. 따라서 이 샘플의 평균은 107.9235 이다. 그런데 연구자는 이 샘플이 나온 모집단은 전혀 모르는 상태이다. 모집단의 평균도 모르고 분산값도 모른다. 그냥 10명을 샘플로 뽑았고, 이들이 p1에서 나왔을 확률을 알아보려고 한다. 이를 알아보기 위해서 스크립트를 돌려 보니, blue와 green라인 밖의 범위가 나올 확률은 0.01222302로 5/100 보다 작고, 1/100 에 가깝다.  
 + 
 +이 probability level이 어느정도나 작아야 이 샘플이 p1에서 나오지 않고 p2에서 나왔다고 판단할 수 있을까? 관습적으로 5/100를 기준으로 해서 이 범위보다 작게 되면 p1의 모집단에서 나온 샘플이 아닌 것으로 판단하게 된다. 이 논리에 따라서, (평균을 107.9235 값을 갖는)이 샘플은, p1에서 나오기 힘든 샘플이라고 판단 된다.  
 + 
 +  * R에는 z-test 펑션이 없다. 현실에서는 전체 모집단의 평균을 알고 있는 경우는 많지만 표준편차까지 알고 있는 경우는 많지 않다. 그래서 많이 쓰이지 않는 편이다.  
 +  * 모집단의 평균과 표준편차를 알고 있다고 하면, 우리는 R에서 z test를 하는 절차는  
 +  * n = n 일 경우의 샘플링분포에서 se 를 구한다 
 +    * se = sigma / sqrt(n) 
 +  * 테스트할 점수의 z score를 구한다. 
 +    * diff = test.score - mean.of.population 
 +    * z.score = diff / se  
 +  * z score 보다 큰 점수나 -z score 보다 작은 점수가 나올 확률를 위의 샘플링 분포에 구한다. 
 +    * p.value = 2 * pnorm(z.score, lower.tail=F) 
 +  * z.score와 p.vallue로 테스트점수가 모집단에서 나왔는지 나올 수 없는지 (나오기 어려운지를) 판단한다.
 {{:r:pasted:20250910-135442.png?400}} {{:r:pasted:20250910-135442.png?400}}
  
Line 868: Line 910:
 <WRAP column half> <WRAP column half>
 <code> <code>
 +> ###################################################
 +> # 다른 케이스 --> t-test 케이스
 > # what if we do not know the sd of  > # what if we do not know the sd of 
 > # population? > # population?
Line 876: Line 920:
 > # of normal distribution (z). > # of normal distribution (z).
 > diff > diff
-[1] 7.95152+[1] 7.923527
 > se.t <- sqrt(var(sample.i.got.p2)/s.size) > se.t <- sqrt(var(sample.i.got.p2)/s.size)
 > t.cal <- diff/se.t > t.cal <- diff/se.t
 > t.cal > t.cal
-[1] 2.417997+[1] 2.914599
 > pt(t.cal, df=s.size-1, lower.tail = F) > pt(t.cal, df=s.size-1, lower.tail = F)
-[1] 0.01936882+[1] 0.008591086
 > pt(t.cal, df=s.size-1, lower.tail = F)*2 > pt(t.cal, df=s.size-1, lower.tail = F)*2
-[1] 0.03873764+[1] 0.01718217
 > prob.t <- pt(t.cal, df=s.size-1, lower.tail = F)*2 > prob.t <- pt(t.cal, df=s.size-1, lower.tail = F)*2
 >  > 
Line 890: Line 934:
 > z.cal > z.cal
          [,1]          [,1]
-[1,] 2.514492+[1,] 2.505639
 > se.z > se.z
          [,1]          [,1]
Line 896: Line 940:
 > prob.z > prob.z
            [,1]            [,1]
-[1,] 0.01192042+[1,] 0.01222302
 >  > 
 > t.cal > t.cal
-[1] 2.417997+[1] 2.914599
 > se.t > se.t
-[1] 3.288473+[1] 2.718566
 > prob.t > prob.t
-[1] 0.03873764+[1] 0.01718217
 >  > 
 > t.test(sample.i.got.p2, mu=mean(p1)) > t.test(sample.i.got.p2, mu=mean(p1))
Line 910: Line 954:
  
 data:  sample.i.got.p2 data:  sample.i.got.p2
-t = 2.418, df = 9, p-value = 0.03874+t = 2.9146, df = 9, p-value = 0.01718
 alternative hypothesis: true mean is not equal to 100 alternative hypothesis: true mean is not equal to 100
 95 percent confidence interval: 95 percent confidence interval:
- 100.5125 115.3906+ 101.7737 114.0733
 sample estimates: sample estimates:
 mean of x  mean of x 
- 107.9515 + 107.9235  
 + 
 + 
 </code> </code>
 </WRAP> </WRAP>
 <WRAP half column> <WRAP half column>
 [[:t-test]] 중에서 2번째 케이스 [[:t-test]] 중에서 2번째 케이스
 +  * 모집단의 평균은 알지만 표준편차 정보는 없는 경우이다. 
 +  * 똑같은 논리로 생각을 해서 모집단의 샘플링분포를 (distribution of sample means) 머리에 두고
 +  * se값을 구한다. 이 때의 se 값은 
 +    * ''se.cal <- sqrt( s^2/n )''값으로 구한다. 
 +    * sigma 대신에 s를 사용한 것에 주목
 +  * z.score에 해당하는 t.score를 구한다. 
 +    * 테스트점수와 모집단 평균의 차이를 구한 후 (''diff = test.score - mean(p1)'')
 +    * se.cal 값으로 나눠준다. 
 +    * ''t.cal <- diff / se.cal''
 +  * t.cal 값 이상, 반대편 점수의 이하가 나올 확률을 구한다. 
 +    * 이 때, 모집단의 표준편차를 사용해서 z.score를 구하지 않았으므로 
 +    * 그리고, 이 probability는 샘플의 크기 n에 영향을 받으므로 n의 크기에 따라서 변화하는 probability distribution을 사용한다. 
 +      * p.value <- pt(t.score, df=n-1, lower.tail=F) * 2
 +  * t.cal과 p.value로 테스트점수가 나올 가능성을 판단하여 가설을 기각하거나 채택한다 (검증한다). 
 </WRAP> </WRAP>
 </WRAP> </WRAP>
  
 <WRAP group> <WRAP group>
 +[[:types of error]]
 +[[:hypothesis testing]]
 +
 <WRAP half column> <WRAP half column>
 <code> <code>
->  +> # 그렇다면 
-> # 그렇다면 모집단에 대한 파라미터 없이+
 > # 두 샘플의 평균과 표준편차만 알고 > # 두 샘플의 평균과 표준편차만 알고
 > # 모집단들의 파라미터는 모를 때는? > # 모집단들의 파라미터는 모를 때는?
Line 939: Line 1001:
 > s1 <- sample(p1, s.size, replace = T) > s1 <- sample(p1, s.size, replace = T)
 > s2 <- sample(p2, s.size, replace = T) > s2 <- sample(p2, s.size, replace = T)
 +> s1
 + [1] 100.46855  97.67673  94.48868  97.51996  94.28300 113.78623  95.44617
 + [8]  92.91108  97.46601 100.66815 113.94452  81.20774  93.77529  93.64388
 +[15] 111.43702 104.34787  98.82403 100.93909 111.08048 110.12275  94.68502
 +[22] 100.14734  80.23712  88.02110 118.89602  99.46530 106.92325  97.13274
 +[29] 106.40479 107.44831
 +> s2
 + [1] 120.28360 105.80657 116.28974 106.63910 108.66632 107.51868 123.66436
 + [8]  96.06523 111.53770 108.78226 115.22628 104.60990 105.07035 114.90022
 +[15] 103.92860 105.06587 113.77416 108.07776 111.92824 112.01554 104.34549
 +[22] 100.51113  84.75464  96.36697 117.09381 116.29935 117.20402 109.61657
 +[29] 106.46650 105.86733
 +> 
 > mean(s1) > mean(s1)
-[1] 98.97605+[1] 100.1133
 > mean(s2) > mean(s2)
-[1] 110.8606+[1] 108.6125
 > sd(s1) > sd(s1)
-[1] 10.29997+[1] 9.132472
 > sd(s2) > sd(s2)
-[1] 12.28259+[1] 7.916677
 >  > 
 > ss(s1) > ss(s1)
-[1] 3076.589+[1] 2418.659
 > ss(s2) > ss(s2)
-[1] 4374.997+[1] 1817.54
 > df1 <- s.size - 1 > df1 <- s.size - 1
 > df2 <- s.size - 1 > df2 <- s.size - 1
Line 959: Line 1034:
 > t.ts.cal <- diff/se.tst > t.ts.cal <- diff/se.tst
 > t.ts.cal > t.ts.cal
-[1] 4.060845+[1] 3.851705
 > pt(t.ts.cal, df=df1+df2, lower.tail = F)*2 > pt(t.ts.cal, df=df1+df2, lower.tail = F)*2
-[1] 0.0001484894+[1] 0.0002954557
 >  > 
 > t.test(s2, s1, var.equal=T) > t.test(s2, s1, var.equal=T)
Line 968: Line 1043:
  
 data:  s2 and s1 data:  s2 and s1
-t = 4.0608, df = 58, p-value = 0.0001485+t = 3.8517, df = 58, p-value = 0.0002955
 alternative hypothesis: true difference in means is not equal to 0 alternative hypothesis: true difference in means is not equal to 0
 95 percent confidence interval: 95 percent confidence interval:
-  6.026262 17.742752+  4.082229 12.916309
 sample estimates: sample estimates:
 mean of x mean of y  mean of x mean of y 
-110.86056  98.97605 + 108.6125  100.1133 
  
->  
 >  > 
 </code> </code>
Line 984: Line 1058:
 </WRAP> </WRAP>
 </WRAP> </WRAP>
-====== Hypothesis e.g ====== 
  
-<WRAP group> 
-위에서 p2의 parameter에 대해서 잘 모른다는 점에 주목하라. 그리고 아래 시나리오를 상상하라.  
-어느 한 모집단의 IQ 평균이 100 이고 표준편차가 10 임을 알고 있다. 확률과통계 교수는 머리가 좋아지는 약을 개발하여 이를 팔아보려고 하고 있다. 이를 위해서 확통교수는 25명의 학생에게 머리가 좋아지는 약을 복용하도록 한 후에 IQ를 측정하였다. 그런데, 그 IQ 평균이 105.45 이다. 이 점수를 가지고 약의 효과가 있는지 검증을 해보력고 한다.  
-<tabbox rs.eg01> 
-<code> 
-m.s <- 105.45 
-m.p <- 100 
-sd.p <- 10 
-v.p <- 100 
-s.sz <- 25 
-s <- rnorm2(s.sz, m.s, sd.p) 
-s 
-se <- sqrt(v.p/s.sz) 
-se 
-diff <- m.s - m.p  
-diff 
-t.cal <- diff / se 
-df.s <- s.sz-1 
-2*pt(t.cal, df=df.s,lower.tail = F) 
-t.test(s, mu=m.p) 
-</code> 
-<tabbox ro.eg01> 
-<code> 
-> m.s <- 105.45 
-> m.p <- 100 
-> sd.p <- 10 
-> v.p <- 100 
-> s.sz <- 25 
-> s <- rnorm2(s.sz, m.s, sd.p) 
-> s 
-           [,1] 
- [1,] 102.84774 
- [2,] 114.05839 
- [3,] 101.95736 
- [4,] 122.88205 
- [5,] 106.27718 
- [6,] 116.96923 
- [7,] 101.78039 
- [8,]  81.37328 
- [9,] 112.28533 
-[10,] 108.26472 
-[11,] 100.31204 
-[12,]  89.84322 
-[13,] 109.29927 
-[14,] 116.19203 
-[15,]  98.35046 
-[16,] 105.84416 
-[17,] 104.34385 
-[18,] 119.76468 
-[19,] 104.44693 
-[20,]  93.43272 
-[21,]  95.42572 
-[22,] 107.99983 
-[23,] 119.88284 
-[24,]  95.36440 
-[25,] 107.05218 
-attr(,"scaled:center") 
-[1] -0.2689983 
-attr(,"scaled:scale") 
-[1] 0.7604765 
-> se <- sqrt(v.p/s.sz) 
-> se 
-[1] 2 
-> diff <- m.s - m.p  
-> diff 
-[1] 5.45 
-> t.cal <- diff / se 
-> df.s <- s.sz-1 
-> 2*pt(t.cal, df=df.s,lower.tail = F) 
-[1] 0.01180885 
-> t.test(s, mu=m.p) 
- 
- One Sample t-test 
- 
-data:  s 
-t = 2.725, df = 24, p-value = 0.01181 
-alternative hypothesis: true mean is not equal to 100 
-95 percent confidence interval: 
- 101.3222 109.5778 
-sample estimates: 
-mean of x  
-   105.45  
- 
->  
-</code> 
-</tabbox> 
-</WRAP> 
  
 ====== R script and output ====== ====== R script and output ======
Line 1325: Line 1311:
 c(hi3, hiz3) c(hi3, hiz3)
  
- +means.10 <- rnorm2(iter, mean(p1), c(se.z)) 
-hist(means, breaks=50, +hist(means.10, breaks=50, 
-     xlim = c(mean(means)-5*sd(means), mean(means)+10*sd(means)), +     xlim = c(mean(p1)-5*se.z, mean(p1)+10*se.z), 
      col = rgb(1, 1, 1, .5))      col = rgb(1, 1, 1, .5))
-abline(v=mean(means), col="black", lwd=3)+abline(v=mean(p1), col="black", lwd=3)
 # abline(v=mean(p2), colo='darkgreen', lwd=3) # abline(v=mean(p2), colo='darkgreen', lwd=3)
-abline(v=c(lo1, hi1, lo2, hi2, lo3, hi3), +abline(v=c(loz1, hiz1, loz2, hiz2, loz3, hiz3), 
        col=c("darkgreen","darkgreen", "blue", "blue", "orange", "orange"),         col=c("darkgreen","darkgreen", "blue", "blue", "orange", "orange"), 
        lwd=2)        lwd=2)
Line 1346: Line 1332:
 m.sample.i.got m.sample.i.got
  
-hist(means, breaks=30,  +hist(means.10, breaks=30,  
-     xlim = c(mean(means)-7*sd(means), mean(means)+10*sd(means)), +     xlim = c(mean(p1)-7*se.z, mean(p1)+10*se.z), 
      col = rgb(1, 1, 1, .5))      col = rgb(1, 1, 1, .5))
-abline(v=mean(means), col="black", lwd=3)+abline(v=mean(p1), col="black", lwd=3)
 abline(v=m.sample.i.got, col='darkgreen', lwd=3) abline(v=m.sample.i.got, col='darkgreen', lwd=3)
  
Line 1356: Line 1342:
 # m.sample.i.got? # m.sample.i.got?
 m.sample.i.got m.sample.i.got
-pnorm(m.sample.i.got, mean(means), sd(means), lower.tail = F) 
 pnorm(m.sample.i.got, mean(p1), se.z, lower.tail = F) pnorm(m.sample.i.got, mean(p1), se.z, lower.tail = F)
  
Line 1364: Line 1349:
 # 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))+tmp <- mean(p1) - (m.sample.i.got - mean(p1)) 
 +tmp
 abline(v=tmp, col='red', lwd=3) abline(v=tmp, col='red', lwd=3)
-2 * pnorm(m.sample.i.got, mean(p1), se.z, lower.tail = F) 
 m.sample.i.got m.sample.i.got
 +# 붉은 선의 왼쪽 부분을 구해서 2를 곱한 값의 범위
 +prob.of.m.sample.i.got <- 2 * pnorm(m.sample.i.got, mean(p1), se.z, lower.tail = F)
 +prob.of.m.sample.i.got
  
-2 * pnorm(m.sample.i.got, mean(p1), se.z, lower.tail = F) 
 (m.sample.i.got - mean(p1))/se.z (m.sample.i.got - mean(p1))/se.z
 z.score <- (m.sample.i.got - mean(p1))/se.z z.score <- (m.sample.i.got - mean(p1))/se.z
 pnorm(z.score, lower.tail = FALSE) pnorm(z.score, lower.tail = FALSE)
 2 * pnorm(z.score, lower.tail = FALSE) 2 * pnorm(z.score, lower.tail = FALSE)
 +
 +# 즉 n = ssize 크기의 샘플을 구했을 때 
 +# 104.7489 점 오른쪽 그리고,
 +# 이에 대응하는 95.25113 지점 왼쪽의 
 +# 평균이 나오는 probability는 
 +# 0.1331684 
  
 ### one more time  ### one more time 
 # this time, with a story # this time, with a story
 +# 한편 우리는 특정한 평균과 표준편차를 갖는
 +# p2를 알고 있다 (110, 10)
 mean(p2) mean(p2)
 sd(p2) sd(p2)
Line 1386: Line 1381:
 m.sample.i.got.p2 m.sample.i.got.p2
  
-hist(means, breaks=30,  +# 위의 샘플 평균이 위에서 얘기했던 샘플평균들의 집합에서  
-     xlim = c(tmp-10*sd(means), m.sample.i.got+10*sd(means)), +# 나올 확률은 무엇일까? 
 +hist(means.10, breaks=30,  
 +     xlim = c(tmp-10*se.z, m.sample.i.got+10*se.z), 
      col = rgb(1, 1, 1, .5))      col = rgb(1, 1, 1, .5))
-abline(v=mean(means), col="black", lwd=3)+abline(v=mean(p1), col="black", lwd=3)
 abline(v=m.sample.i.got.p2, col='blue', lwd=3) abline(v=m.sample.i.got.p2, col='blue', lwd=3)
  
Line 1396: Line 1393:
 # m.sample.i.got? # m.sample.i.got?
 m.sample.i.got.p2 m.sample.i.got.p2
-pnorm(m.sample.i.got.p2, mean(p1), se.z, lower.tail = F)+pnorm(m.sample.i.got.p2, mean(p1), se.z, lower.tail = F)  
 +# 혹은 
 +z.cal <- (m.sample.i.got.p2-mean(p1))/se.z 
 +pnorm(z.cal, lower.tail = F)
  
-# then, what is the probabilty of getting +# then, what is the probability of getting 
 # greater than m.sample.i.got and # greater than m.sample.i.got and
 # less than corresponding value, which is # less than corresponding value, which is
Line 1405: Line 1405:
 abline(v=mean(p1)-(m.sample.i.got.p2-mean(p1)), col='red', lwd=3) abline(v=mean(p1)-(m.sample.i.got.p2-mean(p1)), col='red', lwd=3)
 2 * pnorm(m.sample.i.got.p2, mean(p1), se.z, lower.tail = F) 2 * pnorm(m.sample.i.got.p2, mean(p1), se.z, lower.tail = F)
 +2 * pnorm(z.cal, lower.tail = F)
  
-z.cal <- (m.sample.i.got.p2 - mean(p1)) / se.z +# 위는 샘플평균과 모집단 평균값 간의 차이 diff를 
-se.z +# standard error값으로 즉,  
-z.cal +# sd of sample means (n=s.size) 값으로 나눠 준  
- +# 표준점수 값을 구하고  
-pnorm(z.cal, lower.tail = F) +# 이 점수에 대한 probability값을 구해 본것
-2 * pnorm(z.cal, lower.tail = F)+
  
 diff <- m.sample.i.got.p2 - mean(p1) diff <- m.sample.i.got.p2 - mean(p1)
Line 1417: Line 1417:
 z.cal <- diff / se.z z.cal <- diff / se.z
 z.cal z.cal
-pnorm(z.cal, lower.tail=F)*2 
 prob.z <- pnorm(z.cal, lower.tail=F)*2 prob.z <- pnorm(z.cal, lower.tail=F)*2
 +prob.z
 +# 위에 대한 해석 . . . . . 
 +#
 +#
 +###################################################
  
 +###################################################
 +# 다른 케이스 --> t-test 케이스
 # what if we do not know the sd of  # what if we do not know the sd of 
 # population? # population?
Line 1457: Line 1463:
 s1 <- sample(p1, s.size, replace = T) s1 <- sample(p1, s.size, replace = T)
 s2 <- sample(p2, s.size, replace = T) s2 <- sample(p2, s.size, replace = T)
 +s1
 +s2
 +
 mean(s1) mean(s1)
 mean(s2) mean(s2)
Line 1474: Line 1483:
  
 t.test(s2, s1, var.equal=T) t.test(s2, s1, var.equal=T)
 +
 +
 +comb <- list(s1 = s1, s2 = s2)
 +op <- stack(comb)
 +head(op)
 +colnames(op)[1] <- "values"
 +colnames(op)[2] <- "group"
 +op$group <- factor(op$group)
 +head(op)
 +boxplot(op$values~op$group)
 +
 +# Calculate the mean of 'Value' grouped by 'Category'
 +m.by.group <- aggregate(values ~ group, data = op, FUN = mean, na.rm = TRUE)
 +m.by.group
 +m.s1 <- m.by.group$values[1]
 +m.s2 <- m.by.group$values[2]
 +m.tot <- mean(op$values)
 +m.s1 <- mean(s1)
 +m.s2 <- mean(s2)
 +
 +n.s1 <- length(s1)
 +n.s2 <- length(s2)
 +df.s1 <- n.s1-1
 +df.s2 <- n.s2-1
 +
 +ss.tot <- sum((m.tot-op$values)^2)
 +ss.tot
 +ss.tot <- ss(op$values)
 +ss.tot
 +ss.bet <- ((m.tot-m.s1)^2*n.s1) + 
 +  ((m.tot-m.s2)^2*n.s2)
 +ss.wit <- ss(s1) + ss(s2)
 +ss.tot
 +ss.bet
 +ss.wit
 +ss.bet+ss.wit
 +
 +df.tot <- (n.s1+n.s2)-1
 +df.bet <- nlevels(op$group)-1
 +df.wit <- df.s1+df.s2
 +
 +ms.tot <- ss.tot / df.tot
 +ms.bet <- ss.bet / df.bet
 +ms.wit <- ss.wit / df.wit
 +
 +r.sq <- ss.bet / ss.tot
 +r.sq
 +f.cal <- ms.bet / ms.wit
 +f.cal
 +prob.f <- pf(f.cal, df.bet, df.wit, lower.tail = F)
 +prob.f
 +
 +ss.bet
 +df.bet
 +ms.bet
 +ss.wit
 +df.wit
 +ms.wit
 +f.cal
 +prob.f
 +
 +m.f <- aov(values~group, data=op)
 +summary(m.f)
 +
 +m.t <- t.test(s2, s1, var.equal = T)
 +m.t
 +m.t$statistic^2
 +m.t$p.value
 +
 +m.l <- lm(values~group, data=op)
 +summary(m.l)
 +f.cal
 +df.bet
 +df.wit
 +prob.f
 +r.sq
 +anova(m.l)
 +ss.bet
 +ss.wit
 +ms.bet
 +ms.wit
 +f.cal
 +prob.f
 </code> </code>
  
Line 1579: Line 1671:
 >  > 
 > mean(za.p1) > mean(za.p1)
-[1] 2.884512e-17+[1] -2.100958e-17
 > mean(za.p1)==mean(zb.p1) > mean(za.p1)==mean(zb.p1)
 [1] TRUE [1] TRUE
Line 1711: Line 1803:
 > s1 <- sample(p1, s.size, replace = T) > s1 <- sample(p1, s.size, replace = T)
 > mean(s1) > mean(s1)
-[1] 100.7648+[1] 100.7771
 > means <- append(means, mean(s1)) > means <- append(means, mean(s1))
 > means <- append(means, mean(sample(p1, s.size, replace = T))) > means <- append(means, mean(sample(p1, s.size, replace = T)))
Line 1718: Line 1810:
 > means <- append(means, mean(sample(p1, s.size, replace = T))) > means <- append(means, mean(sample(p1, s.size, replace = T)))
 > means > means
-[1] 100.7648 105.3139 107.5766 102.4201 101.6467+[1] 100.77711  99.34268 101.12569  99.03624  98.73117
 > hist(means) > hist(means)
 >  > 
Line 1734: Line 1826:
 [1] 100000 [1] 100000
 > mean(means) > mean(means)
-[1] 100.0025+[1] 99.99036
 > sd(means) > sd(means)
-[1] 3.167479+[1] 3.156405
 > var(means) > var(means)
-[1] 10.03292+[1] 9.962894
 > se.s <- sd(means) > se.s <- sd(means)
 >  > 
Line 1760: Line 1852:
 > # meanwhile . . . . > # meanwhile . . . .
 > se.s > se.s
-[1] 3.167479+[1] 3.156405
 > sd(means) > sd(means)
-[1] 3.167479+[1] 3.156405
 > var(means) > var(means)
-[1] 10.03292+[1] 9.962894
 >  > 
 > se.z <- sqrt(var(p1)/s.size) > se.z <- sqrt(var(p1)/s.size)
Line 1789: Line 1881:
 >  > 
 > mean(means) > mean(means)
-[1] 100.0025+[1] 99.99036
 > var(means)  > var(means) 
-[1] 10.03292+[1] 9.962894
 > sd(means) > sd(means)
-[1] 3.167479+[1] 3.156405
 > # remember we started talking sample size 10 > # remember we started talking sample size 10
 > mean(p1) > mean(p1)
Line 1813: Line 1905:
 >  > 
 > c(lo1, loz1) > c(lo1, loz1)
-[1] 96.83502 96.83772+[1] 96.83395 96.83772
 > c(lo2, loz2) > c(lo2, loz2)
-[1] 93.66754 93.67544+[1] 93.67755 93.67544
 > c(lo3, loz3) > c(lo3, loz3)
-[1] 90.50006 90.51317+[1] 90.52114 90.51317
 >  > 
 > c(hi1, hiz1) > c(hi1, hiz1)
-[1] 103.1700 103.1623+[1] 103.1468 103.1623
 > c(hi2, hiz2) > c(hi2, hiz2)
-[1] 106.3375 106.3246+[1] 106.3032 106.3246
 > c(hi3, hiz3) > c(hi3, hiz3)
-[1] 109.5049 109.4868+[1] 109.4596 109.4868
 >  > 
->  +> means.10 <- rnorm2(iter, mean(p1), c(se.z)) 
-> hist(means, breaks=50, +> hist(means.10, breaks=50, 
-+      xlim = c(mean(means)-5*sd(means), mean(means)+10*sd(means)), ++      xlim = c(mean(p1)-5*se.z, mean(p1)+10*se.z), 
 +      col = rgb(1, 1, 1, .5)) +      col = rgb(1, 1, 1, .5))
-> abline(v=mean(means), col="black", lwd=3)+> abline(v=mean(p1), col="black", lwd=3)
 > # abline(v=mean(p2), colo='darkgreen', lwd=3) > # abline(v=mean(p2), colo='darkgreen', lwd=3)
-> abline(v=c(lo1, hi1, lo2, hi2, lo3, hi3), +> abline(v=c(loz1, hiz1, loz2, hiz2, loz3, hiz3), 
 +        col=c("darkgreen","darkgreen", "blue", "blue", "orange", "orange"),  +        col=c("darkgreen","darkgreen", "blue", "blue", "orange", "orange"), 
 +        lwd=2) +        lwd=2)
Line 1841: Line 1933:
 [1]  94 106 [1]  94 106
 > round(c(lo3, hi3)) > round(c(lo3, hi3))
-[1]  91 110+[1]  91 109
 >  > 
 > round(c(loz1, hiz1)) > round(c(loz1, hiz1))
Line 1852: Line 1944:
 > m.sample.i.got <- mean(means)+ 1.5*sd(means) > m.sample.i.got <- mean(means)+ 1.5*sd(means)
 > m.sample.i.got > m.sample.i.got
-[1] 104.7537+[1] 104.725
 >  > 
-> hist(means, breaks=30,  +> hist(means.10, breaks=30,  
-+      xlim = c(mean(means)-7*sd(means), mean(means)+10*sd(means)), ++      xlim = c(mean(p1)-7*se.z, mean(p1)+10*se.z), 
 +      col = rgb(1, 1, 1, .5)) +      col = rgb(1, 1, 1, .5))
-> abline(v=mean(means), col="black", lwd=3)+> abline(v=mean(p1), col="black", lwd=3)
 > abline(v=m.sample.i.got, col='darkgreen', lwd=3) > abline(v=m.sample.i.got, col='darkgreen', lwd=3)
 >  > 
Line 1864: Line 1956:
 > # m.sample.i.got? > # m.sample.i.got?
 > m.sample.i.got > m.sample.i.got
-[1] 104.7537 +[1] 104.725
-> pnorm(m.sample.i.got, mean(means), sd(means), lower.tail = F) +
-[1] 0.0668072+
 > pnorm(m.sample.i.got, mean(p1), se.z, lower.tail = F) > pnorm(m.sample.i.got, mean(p1), se.z, lower.tail = F)
-[1] 0.06638638+[1] 0.06756624
 >  > 
 > # then, what is the probabilty of getting  > # then, what is the probabilty of getting 
Line 1875: Line 1965:
 > # 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))+> tmp <- mean(p1) - (m.sample.i.got - mean(p1)) 
 +> tmp 
 +[1] 95.27504
 > abline(v=tmp, col='red', lwd=3) > abline(v=tmp, col='red', lwd=3)
-> 2 * pnorm(m.sample.i.got, mean(p1), se.z, lower.tail = F) 
-[1] 0.1327728 
 > m.sample.i.got > m.sample.i.got
-[1] 104.7537+[1] 104.725 
 +> # 붉은 선의 왼쪽 부분을 구해서 2를 곱한 값의 범위 
 +> prob.of.m.sample.i.got <- 2 * pnorm(m.sample.i.got, mean(p1), se.z, lower.tail = F) 
 +> prob.of.m.sample.i.got 
 +[1] 0.1351325
 >  > 
-> 2 * pnorm(m.sample.i.got, mean(p1), se.z, lower.tail = F) 
-[1] 0.1327728 
 > (m.sample.i.got - mean(p1))/se.z > (m.sample.i.got - mean(p1))/se.z
          [,1]          [,1]
-[1,] 1.503257+[1,] 1.494165
 > z.score <- (m.sample.i.got - mean(p1))/se.z > z.score <- (m.sample.i.got - mean(p1))/se.z
 > pnorm(z.score, lower.tail = FALSE) > pnorm(z.score, lower.tail = FALSE)
            [,1]            [,1]
-[1,] 0.06638638+[1,] 0.06756624
 > 2 * pnorm(z.score, lower.tail = FALSE) > 2 * pnorm(z.score, lower.tail = FALSE)
           [,1]           [,1]
-[1,] 0.1327728+[1,] 0.1351325 
 +>  
 +> # 즉 n = ssize 크기의 샘플을 구했을 때  
 +> # 104.7489 점 오른쪽 그리고, 
 +> # 이에 대응하는 95.25113 지점 왼쪽의  
 +> # 평균이 나오는 probability는  
 +> # 0.1331684 
 >  > 
 > ### one more time  > ### one more time 
 > # this time, with a story > # this time, with a story
 +> # 한편 우리는 특정한 평균과 표준편차를 갖는
 +> # p2를 알고 있다 (110, 10)
 > mean(p2) > mean(p2)
 [1] 110 [1] 110
Line 1904: Line 2004:
 > sample.i.got.p2 <- sample(p2, s.size) > sample.i.got.p2 <- sample(p2, s.size)
 > sample.i.got.p2 > sample.i.got.p2
- [1] 109.64312 117.28119 106.40247 113.29156 124.23996 117.81039  91.88752 + [1] 106.63828 107.09795 108.08661 119.09268 103.38813 118.30613 111.49542 
- [8] 102.20386  95.59136 101.16378+ [8]  94.25277 115.38045  95.49686
 >  > 
 > m.sample.i.got.p2 <- mean(sample.i.got.p2) > m.sample.i.got.p2 <- mean(sample.i.got.p2)
 > m.sample.i.got.p2 > m.sample.i.got.p2
-[1] 107.9515+[1] 107.9235
 >  > 
-> hist(means, breaks=30,  +> # 위의 샘플 평균이 위에서 얘기했던 샘플평균들의 집합에서  
-+      xlim = c(tmp-10*sd(means), m.sample.i.got+10*sd(means)), +> # 나올 확률은 무엇일까? 
 +> hist(means.10, breaks=30,  
 ++      xlim = c(tmp-10*se.z, m.sample.i.got+10*se.z), 
 +      col = rgb(1, 1, 1, .5)) +      col = rgb(1, 1, 1, .5))
-> abline(v=mean(means), col="black", lwd=3)+> abline(v=mean(p1), col="black", lwd=3)
 > abline(v=m.sample.i.got.p2, col='blue', lwd=3) > abline(v=m.sample.i.got.p2, col='blue', lwd=3)
 >  > 
Line 1921: Line 2023:
 > # m.sample.i.got? > # m.sample.i.got?
 > m.sample.i.got.p2 > m.sample.i.got.p2
-[1] 107.9515 +[1] 107.9235 
-> pnorm(m.sample.i.got.p2, mean(p1), se.z, lower.tail = F) +> pnorm(m.sample.i.got.p2, mean(p1), se.z, lower.tail = F)  
-[1] 0.005960209+[1] 0.006111512 
 +> # 혹은 
 +> z.cal <- (m.sample.i.got.p2-mean(p1))/se.z 
 +> pnorm(z.cal, lower.tail = F) 
 +            [,1] 
 +[1,] 0.006111512
 >  > 
-> # then, what is the probabilty of getting +> # then, what is the probability of getting 
 > # greater than m.sample.i.got and > # greater than m.sample.i.got and
 > # less than corresponding value, which is > # less than corresponding value, which is
Line 1932: Line 2039:
 > abline(v=mean(p1)-(m.sample.i.got.p2-mean(p1)), col='red', lwd=3) > abline(v=mean(p1)-(m.sample.i.got.p2-mean(p1)), col='red', lwd=3)
 > 2 * pnorm(m.sample.i.got.p2, mean(p1), se.z, lower.tail = F) > 2 * pnorm(m.sample.i.got.p2, mean(p1), se.z, lower.tail = F)
-[1] 0.01192042 +[1] 0.01222302
->  +
-> z.cal <- (m.sample.i.got.p2 - mean(p1)) / se.z +
-> se.z +
-         [,1] +
-[1,] 3.162278 +
-> z.cal +
-         [,1] +
-[1,] 2.514492 +
->  +
-> pnorm(z.cal, lower.tail = F) +
-            [,1] +
-[1,] 0.005960209+
 > 2 * pnorm(z.cal, lower.tail = F) > 2 * pnorm(z.cal, lower.tail = F)
            [,1]            [,1]
-[1,] 0.01192042+[1,] 0.01222302 
 +>  
 +> # 위는 샘플평균과 모집단 평균값 간의 차이 diff를 
 +> # standard error값으로 즉,  
 +> # sd of sample means (n=s.size) 값으로 나눠 준  
 +> # 표준점수 값을 구하고  
 +> # 이 점수에 대한 probability값을 구해 본것
 >  > 
 > diff <- m.sample.i.got.p2 - mean(p1) > diff <- m.sample.i.got.p2 - mean(p1)
Line 1954: Line 2055:
 > z.cal > z.cal
          [,1]          [,1]
-[1,] 2.514492 +[1,] 2.505639 
-> pnorm(z.cal, lower.tail=F)*2+> prob.z <- pnorm(z.cal, lower.tail=F)*2 
 +> prob.z
            [,1]            [,1]
-[1,] 0.01192042 +[1,] 0.01222302 
-> prob.z <- pnorm(z.cal, lower.tail=F)*2+> # 위에 대한 해석 . . . . .  
 +> # 
 +> # 
 +> ###################################################
 >  > 
 +> ###################################################
 +> # 다른 케이스 --> t-test 케이스
 > # what if we do not know the sd of  > # what if we do not know the sd of 
 > # population? > # population?
Line 1968: Line 2075:
 > # of normal distribution (z). > # of normal distribution (z).
 > diff > diff
-[1] 7.95152+[1] 7.923527
 > se.t <- sqrt(var(sample.i.got.p2)/s.size) > se.t <- sqrt(var(sample.i.got.p2)/s.size)
 > t.cal <- diff/se.t > t.cal <- diff/se.t
 > t.cal > t.cal
-[1] 2.417997+[1] 2.914599
 > pt(t.cal, df=s.size-1, lower.tail = F) > pt(t.cal, df=s.size-1, lower.tail = F)
-[1] 0.01936882+[1] 0.008591086
 > pt(t.cal, df=s.size-1, lower.tail = F)*2 > pt(t.cal, df=s.size-1, lower.tail = F)*2
-[1] 0.03873764+[1] 0.01718217
 > prob.t <- pt(t.cal, df=s.size-1, lower.tail = F)*2 > prob.t <- pt(t.cal, df=s.size-1, lower.tail = F)*2
 >  > 
Line 1982: Line 2089:
 > z.cal > z.cal
          [,1]          [,1]
-[1,] 2.514492+[1,] 2.505639
 > se.z > se.z
          [,1]          [,1]
Line 1988: Line 2095:
 > prob.z > prob.z
            [,1]            [,1]
-[1,] 0.01192042+[1,] 0.01222302
 >  > 
 > t.cal > t.cal
-[1] 2.417997+[1] 2.914599
 > se.t > se.t
-[1] 3.288473+[1] 2.718566
 > prob.t > prob.t
-[1] 0.03873764+[1] 0.01718217
 >  > 
 > t.test(sample.i.got.p2, mu=mean(p1)) > t.test(sample.i.got.p2, mu=mean(p1))
Line 2002: Line 2109:
  
 data:  sample.i.got.p2 data:  sample.i.got.p2
-t = 2.418, df = 9, p-value = 0.03874+t = 2.9146, df = 9, p-value = 0.01718
 alternative hypothesis: true mean is not equal to 100 alternative hypothesis: true mean is not equal to 100
 95 percent confidence interval: 95 percent confidence interval:
- 100.5125 115.3906+ 101.7737 114.0733
 sample estimates: sample estimates:
 mean of x  mean of x 
- 107.9515 + 107.9235 
  
 >  > 
Line 2022: Line 2129:
 > s1 <- sample(p1, s.size, replace = T) > s1 <- sample(p1, s.size, replace = T)
 > s2 <- sample(p2, s.size, replace = T) > s2 <- sample(p2, s.size, replace = T)
 +> s1
 + [1] 100.46855  97.67673  94.48868  97.51996  94.28300 113.78623  95.44617
 + [8]  92.91108  97.46601 100.66815 113.94452  81.20774  93.77529  93.64388
 +[15] 111.43702 104.34787  98.82403 100.93909 111.08048 110.12275  94.68502
 +[22] 100.14734  80.23712  88.02110 118.89602  99.46530 106.92325  97.13274
 +[29] 106.40479 107.44831
 +> s2
 + [1] 120.28360 105.80657 116.28974 106.63910 108.66632 107.51868 123.66436
 + [8]  96.06523 111.53770 108.78226 115.22628 104.60990 105.07035 114.90022
 +[15] 103.92860 105.06587 113.77416 108.07776 111.92824 112.01554 104.34549
 +[22] 100.51113  84.75464  96.36697 117.09381 116.29935 117.20402 109.61657
 +[29] 106.46650 105.86733
 +> 
 > mean(s1) > mean(s1)
-[1] 98.97605+[1] 100.1133
 > mean(s2) > mean(s2)
-[1] 110.8606+[1] 108.6125
 > sd(s1) > sd(s1)
-[1] 10.29997+[1] 9.132472
 > sd(s2) > sd(s2)
-[1] 12.28259+[1] 7.916677
 >  > 
 > ss(s1) > ss(s1)
-[1] 3076.589+[1] 2418.659
 > ss(s2) > ss(s2)
-[1] 4374.997+[1] 1817.54
 > df1 <- s.size - 1 > df1 <- s.size - 1
 > df2 <- s.size - 1 > df2 <- s.size - 1
Line 2042: Line 2162:
 > t.ts.cal <- diff/se.tst > t.ts.cal <- diff/se.tst
 > t.ts.cal > t.ts.cal
-[1] 4.060845+[1] 3.851705
 > pt(t.ts.cal, df=df1+df2, lower.tail = F)*2 > pt(t.ts.cal, df=df1+df2, lower.tail = F)*2
-[1] 0.0001484894+[1] 0.0002954557
 >  > 
 > t.test(s2, s1, var.equal=T) > t.test(s2, s1, var.equal=T)
Line 2051: Line 2171:
  
 data:  s2 and s1 data:  s2 and s1
-t = 4.0608, df = 58, p-value = 0.0001485+t = 3.8517, df = 58, p-value = 0.0002955
 alternative hypothesis: true difference in means is not equal to 0 alternative hypothesis: true difference in means is not equal to 0
 95 percent confidence interval: 95 percent confidence interval:
-  6.026262 17.742752+  4.082229 12.916309
 sample estimates: sample estimates:
 mean of x mean of y  mean of x mean of y 
-110.86056  98.97605 + 108.6125  100.1133 
  
 +> 
 +> 
 +> comb <- list(s1 = s1, s2 = s2)
 +> op <- stack(comb)
 +> head(op)
 +     values ind
 +1 100.46855  s1
 +2  97.67673  s1
 +3  94.48868  s1
 +4  97.51996  s1
 +5  94.28300  s1
 +6 113.78623  s1
 +> colnames(op)[1] <- "values"
 +> colnames(op)[2] <- "group"
 +> op$group <- factor(op$group)
 +> head(op)
 +     values group
 +1 100.46855    s1
 +2  97.67673    s1
 +3  94.48868    s1
 +4  97.51996    s1
 +5  94.28300    s1
 +6 113.78623    s1
 +> boxplot(op$values~op$group)
 +> 
 +> # Calculate the mean of 'Value' grouped by 'Category'
 +> m.by.group <- aggregate(values ~ group, data = op, FUN = mean, na.rm = TRUE)
 +> m.by.group
 +  group   values
 +1    s1 100.1133
 +2    s2 108.6125
 +> m.s1 <- m.by.group$values[1]
 +> m.s2 <- m.by.group$values[2]
 +> m.tot <- mean(op$values)
 +> m.s1 <- mean(s1)
 +> m.s2 <- mean(s2)
 +> 
 +> n.s1 <- length(s1)
 +> n.s2 <- length(s2)
 +> df.s1 <- n.s1-1
 +> df.s2 <- n.s2-1
 +> 
 +> ss.tot <- sum((m.tot-op$values)^2)
 +> ss.tot
 +[1] 5319.762
 +> ss.tot <- ss(op$values)
 +> ss.tot
 +[1] 5319.762
 +> ss.bet <- ((m.tot-m.s1)^2*n.s1) + 
 ++   ((m.tot-m.s2)^2*n.s2)
 +> ss.wit <- ss(s1) + ss(s2)
 +> ss.tot
 +[1] 5319.762
 +> ss.bet
 +[1] 1083.564
 +> ss.wit
 +[1] 4236.199
 +> ss.bet+ss.wit
 +[1] 5319.762
 +> 
 +> df.tot <- (n.s1+n.s2)-1
 +> df.bet <- nlevels(op$group)-1
 +> df.wit <- df.s1+df.s2
 +> 
 +> ms.tot <- ss.tot / df.tot
 +> ms.bet <- ss.bet / df.bet
 +> ms.wit <- ss.wit / df.wit
 +> 
 +> r.sq <- ss.bet / ss.tot
 +> r.sq
 +[1] 0.2036865
 +> f.cal <- ms.bet / ms.wit
 +> f.cal
 +[1] 14.83563
 +> prob.f <- pf(f.cal, df.bet, df.wit, lower.tail = F)
 +> prob.f
 +[1] 0.0002954557
 +> 
 +> ss.bet
 +[1] 1083.564
 +> df.bet
 +[1] 1
 +> ms.bet
 +[1] 1083.564
 +> ss.wit
 +[1] 4236.199
 +> df.wit
 +[1] 58
 +> ms.wit
 +[1] 73.03791
 +> f.cal
 +[1] 14.83563
 +> prob.f
 +[1] 0.0002954557
 +> 
 +> m.f <- aov(values~group, data=op)
 +> summary(m.f)
 +            Df Sum Sq Mean Sq F value   Pr(>F)    
 +group        1   1084    1084   14.84 0.000295 ***
 +Residuals   58   4236      73                     
 +---
 +Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
 +> 
 +> m.t <- t.test(s2, s1, var.equal = T)
 +> m.t
 +
 + Two Sample t-test
 +
 +data:  s2 and s1
 +t = 3.8517, df = 58, p-value = 0.0002955
 +alternative hypothesis: true difference in means is not equal to 0
 +95 percent confidence interval:
 +  4.082229 12.916309
 +sample estimates:
 +mean of x mean of y 
 + 108.6125  100.1133 
 +
 +> m.t$statistic^2
 +       t 
 +14.83563 
 +> m.t$p.value
 +[1] 0.0002954557
 +> 
 +> m.l <- lm(values~group, data=op)
 +> summary(m.l)
 +
 +Call:
 +lm(formula = values ~ group, data = op)
 +
 +Residuals:
 +     Min       1Q   Median       3Q      Max 
 +-23.8579  -4.3671  -0.5914   6.3721  18.7827 
 +
 +Coefficients:
 +            Estimate Std. Error t value Pr(>|t|)    
 +(Intercept)  100.113      1.560  64.162  < 2e-16 ***
 +groups2        8.499      2.207   3.852 0.000295 ***
 +---
 +Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
 +
 +Residual standard error: 8.546 on 58 degrees of freedom
 +Multiple R-squared:  0.2037, Adjusted R-squared:   0.19 
 +F-statistic: 14.84 on 1 and 58 DF,  p-value: 0.0002955
 +
 +> f.cal
 +[1] 14.83563
 +> df.bet
 +[1] 1
 +> df.wit
 +[1] 58
 +> prob.f
 +[1] 0.0002954557
 +> r.sq
 +[1] 0.2036865
 +> anova(m.l)
 +Analysis of Variance Table
 +
 +Response: values
 +          Df Sum Sq Mean Sq F value    Pr(>F)    
 +group      1 1083.6 1083.56  14.836 0.0002955 ***
 +Residuals 58 4236.2   73.04                      
 +---
 +Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
 +> ss.bet
 +[1] 1083.564
 +> ss.wit
 +[1] 4236.199
 +> ms.bet
 +[1] 1083.564
 +> ms.wit
 +[1] 73.03791
 +> f.cal
 +[1] 14.83563
 +> prob.f
 +[1] 0.0002954557
 >  > 
 </code> </code>
 </tabbox> </tabbox>
  
 +====== exercise.1 ======
 +위에서 p2의 parameter에 대해서 잘 모른다는 점에 주목하라. 그리고 아래 시나리오를 상상하라. 
 +
 +  - 어느 한 모집단의 IQ 평균이 100 이고 표준편차가 10 임을 알고 있다. 확률과통계 교수는 머리가 좋아지는 약을 개발하여 이를 팔아보려고 하고 있다. 이를 위해서 확통교수는 25명의 학생에게 머리가 좋아지는 약을 복용하도록 한 후에 IQ를 측정하였다. 그런데, 그 IQ 평균이 106.45 이다. 이 점수를 가지고 약의 효과가 있는지 검증을 해보력고 한다. 
 +  - 똑 같은 경우이지만 모집단의 평균을 100으로 추정하고 있지 표준편차는 모르는 상태이다. 그런데, 25명에게 약을 복용시키고 IQ를 측정하니 점수가 105.50이 나왔고, 표준편차는 9.4였다. 이 점수로 약의 효과가 있었는지 검증을 하려 한다. 
r/sampling_distribution.1774088927.txt.gz · Last modified: by hkimscil
Donate Powered by PHP Valid HTML5 Valid CSS Driven by DokuWiki