Differences

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

--- r:sampling_distribution [2025/09/10 13:53] – [Last one . . . Important] hkimscil
+++ r:sampling_distribution [2025/09/10 20:41] (current) – [PS1. week02] hkimscil
@@ Line 1: / Line 1: @@
 ====== PS1. week02 ======
 <code>
 rm(list=ls())
@@ Line 32: / Line 31: @@
 var(p1)
-hist(p1, breaks=100, col=rgb(1,1,1,1))
+hist(p1, breaks=50, col = rgb(1, 1, 1, 0.5),
+     main = "histogram of p1 and p2",)
+abline(v=mean(p1), col="black", lwd=3)
+hist(p2, add=T, breaks=50, col=rgb(1,1,.5,.5))
+abline(v=mean(p2), col="red", lwd=3)
+hist(p1, breaks=50, col=rgb(0,.5,.5,.5))
 abline(v=mean(p1),lwd=2)
 abline(v=mean(p1)-sd(p1), lwd=2)
@@ Line 90: / Line 97: @@
 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, -1.8, 1.8), col='red')
 -(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),lwd=2)
 abline(v=mean(p1)-sd(p1), lwd=2)
@@ Line 116: / Line 123: @@
 c(a, b)
 c(-1, 1)
+# note that
+.32/2
+pnorm(-1)
+qnorm(.32/2)
+qnorm(pnorm(-1))
 # 95%
@@ Line 122: / Line 134: @@
 c(c, d)
 c(-2,2)
 # 99%
 e <- qnorm(.01/2)
@@ Line 127: / Line 140: @@
 c(e,f)
 c(-3,3)
 pnorm(b)-pnorm(a)
@@ Line 140: / Line 154: @@
 ################################
-hist(p1, breaks=50, col = rgb(1, 0, 0, 0.5),
-     main = "histogram of p1 and p2",)
-abline(v=mean(p1), col="black", lwd=3)
-hist(p2, add=T, breaks=50, col=rgb(0,0,1,.5))
-abline(v=mean(p2), col="violet", lwd=3)
 s.size <- 10
@@ Line 170: / Line 178: @@
 se.s <- sd(means)
-hist(means, breaks=100, col=rgb(.1, 0, 0, .5))
+hist(means, breaks=50,
-abline(v=mean(means), col="red", lwd=2)
+     xlim = c(mean(means)-5*sd(means), mean(means)+10*sd(means)),
+     col=rgb(1, 1, 1, .5))
+abline(v=mean(means), col="black", lwd=3)
 # now we want to get sd of this distribution
 lo1 <- mean(means)-se.s
@@ Line 180: / Line 189: @@
 lo3 <- mean(means)-3*se.s
 hi3 <- mean(means)+3*se.s
-hist(means,
-     xlim = c(mean(means)-5*sd(means), mean(means)+10*sd(means)),
-     col = rgb(1, 0, 0, .5))
 abline(v=mean(means), col="black", lwd=2)
-# abline(v=mean(p2), colo='darkgreen', lwd=3)
+# abline(v=mean(p2), colo='darkgreen', lwd=2)
 abline(v=c(lo1, hi1, lo2, hi2, lo3, hi3),
-       col=c("green","green", "blue", "blue", "orange", "orange"),
+       col=c("red","red", "blue", "blue", "orange", "orange"),
        lwd=2)
@@ Line 198: / Line 203: @@
 # sd of sample means (sd(means))
-# is sqrt(var(s1)/s.size) or
-# sd(s1) / sqrt(s.size)
 # = se.s
 # when iter value goes to
-# unlimited value:
+# infinite value:
 # mean(means) = mean(p1)
 # and
 # sd(means) = sd(p1) / sqrt(s.size)
-# that is, sd(means) = se.z
+# that is, se.s = se.z
 # This is called CLT (Central Limit Theorem)
+# see http://commres.net/wiki/cetral_limit_theorem
 mean(means)
 mean(p1)
 sd(means)
 var(p1)
+# remember we started talking sample size 10
 sqrt(var(p1)/s.size)
 se.z
@@ Line 237: / Line 243: @@
-hist(means,
+hist(means, breaks=50,
      xlim = c(mean(means)-5*sd(means), mean(means)+10*sd(means)),
-     col = rgb(1, 0, 0, .5))
+     col = rgb(1, 1, 1, .5))
-abline(v=mean(means), col="black", lwd=2)
+abline(v=mean(means), col="black", lwd=3)
 # abline(v=mean(p2), colo='darkgreen', lwd=3)
 abline(v=c(lo1, hi1, lo2, hi2, lo3, hi3),
-       col=c("green","green", "blue", "blue", "orange", "orange"),
+       col=c("darkgreen","darkgreen", "blue", "blue", "orange", "orange"),
        lwd=2)
@@ Line 257: / Line 263: @@
 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, 0, 0, .5))
+     col = rgb(1, 1, 1, .5))
 abline(v=mean(means), col="black", lwd=3)
 abline(v=m.sample.i.got, col='darkgreen', lwd=3)
@@ Line 276: / Line 282: @@
 # (green line)
 tmp <- mean(means) - (m.sample.i.got - mean(means))
-abline(v=tmp, col='green', lwd=3)
+abline(v=tmp, col='red', lwd=3)
 * pnorm(m.sample.i.got, mean(p1), sd(means), lower.tail = F)
 m.sample.i.got
 ### one more time
+# this time, with a story
 mean(p2)
 sd(p2)
@@ Line 287: / Line 294: @@
 m.sample.i.got
-hist(means,
+tmp <- mean(means) - (m.sample.i.got-mean(means))
-     xlim = c(mean(means)-15*sd(means), mean(means)+15*sd(means)),
+tmp
-     col = rgb(1, 0, 0, .5))
-abline(v=mean(means), col="black", lwd=2)
+hist(means, breaks=30,
-abline(v=m.sample.i.got, col='darkgreen', lwd=2)
+     xlim = c(tmp-4*sd(means), m.sample.i.got+4*sd(means)),
+     col = rgb(1, 1, 1, .5))
+abline(v=mean(means), col="black", lwd=3)
+abline(v=m.sample.i.got, col='blue', lwd=3)
 # what is the probablity of getting
@@ Line 304: / Line 314: @@
 # mean(means) - m.sample.i.got - mean(means)
 # (green line)
-tmp <- mean(means) - (m.sample.i.got - mean(means))
+abline(v=tmp, col='red', lwd=3)
-abline(v=tmp, col='green', lwd=2)
 * pnorm(m.sample.i.got, mean(p1), sd(means), lower.tail = F)
 </code>
 ====== output ======
@@ Line 644: / Line 649: @@
 </WRAP>
 <WRAP column half>
-....
+qnorm
   * qnorm는 pnorm의 반대값을 구하는 명령어
   * 히스토그램에서 검정 색 부분의 바깥 쪽 부분은 32%이고 왼 쪽의 것은 이것의 반인 16% 이다.
@@ Line 1064: / Line 1069: @@
 이전 처럼 10명을 샘플로 뽑는데, p2에서 뽑는다. 따라서 이 샘플의 평균은 114.155 이다. 그런데 연구자는 이 샘플이 나온 모집단은 전혀 모르는 상태이다. 모집단의 평균도 모르고 분산값도 모른다. 그냥 10명을 샘플로 뽑았고, 이들이 p1에서 나왔을 확률을 알아보려고 한다. 이를 알아보기 위해서 스크립트를 돌려 보니, blue와 green라인 밖의 범위가 나올 확률은 7.641611e-06로 0에 가깝다. 따라서, 이 샘플은 10명은 p1에 속한 원소가 아닌 다른 성격을 가진 원소이다.
+위에서 p2의 parameter에 대해서 잘 모른다는 점에 주목하라. 그리고 아래 시나리오를 상상하라.
+어느 한 모집단의 IQ 평균이 100 이고 표준편차가 10 임을 알고 있다. 확률과통계 교수는 머리가 좋아지는 약을 개발하여 이를 팔아보려고 하고 있다. 이를 위해서 확통교수는 10을 뽑아서 머리가 좋아지는 약을 복용하도록 한 후에 IQ를 측정하였다. 그런데, 그 IQ 평균이 114.155 이다. 이 점수를 가지고 약의 효과가 있는지 검증을 해보력고 한다.
 </WRAP>
 </WRAP>
-{{:r:pasted:20250910-134658.png}}
+{{:r:pasted:20250910-135442.png}}