Differences

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

--- r:sampling_distribution [2025/09/09 18:04] – hkimscil
+++ r:sampling_distribution [2025/09/10 08:43] (current) – [qnorm] hkimscil
@@ Line 14: / Line 14: @@
   sum((x-mean(x))^2)
 }
 ################################
+N.p <- 1000000
+m.p <- 100
+sd.p <- 10
+p1 <- rnorm2(N.p, m.p, sd.p)
+mean(p1)
+sd(p1)
+p2 <- rnorm2(N.p, m.p+20, sd.p)
+mean(p2)
+sd(p2)
+m.p1 <- mean(p1)
+sd.p1 <- sd(p1)
+var(p1)
+hist(p1, breaks=100, col=rgb(1,1,1,1))
+abline(v=mean(p1),lwd=2)
+abline(v=mean(p1)-sd(p1), lwd=2)
+abline(v=mean(p1)+sd(p1), lwd=2)
+abline(v=c(m.p1-2*sd.p1, m.p1+2*sd.p1), lwd=2, col='red')
+abline(v=c(m.p1-3*sd.p1, m.p1+3*sd.p1), lwd=2, col='green')
+# area bet black = 68%
+# between red = 95%
+# between green = 99%
+pnorm(m.p1+sd.p1, m.p1, sd.p1)
+pnorm(m.p1-sd.p1, m.p1, sd.p1)
+pnorm(m.p1+sd.p1, m.p1, sd.p1) -
+  pnorm(m.p1-sd.p1, m.p1, sd.p1)
+pnorm(m.p1+2*sd.p1, m.p1, sd.p1) -
+  pnorm(m.p1-2*sd.p1, m.p1, sd.p1)
+pnorm(m.p1+3*sd.p1, m.p1, sd.p1) -
+  pnorm(m.p1-3*sd.p1, m.p1, sd.p1)
+m.p1
+sd.p1
+(m.p1-m.p1)/sd.p1
+((m.p1-sd.p1) - m.p1) / sd.p1
+(120-100)/10
+pnorm(1)-pnorm(-1)
+pnorm(2)-pnorm(-2)
+pnorm(3)-pnorm(3)
+-pnorm(-2)*2
+pnorm(2)-pnorm(-2)
+pnorm(120, 100, 10)
+pnorm(2)-pnorm(-2)
+zscore <- (120-100)/10
+pnorm(zscore)-pnorm(-zscore)
+zscore
 # reminder.
 pnorm(-1)
-pnorm(1, lower.tail = F)
+pnorm(1, 0, 1, lower.tail = F)
--(pnorm(-1) + pnorm(1, lower.tail = F))
+pnorm(110, 100, 10, lower.tail = F)
+zscore <- (110-100)/10
+pnorm(zscore, lower.tail = F)
+pnorm(118, 100, 10, lower.tail = F)
+pnorm(18/10, lower.tail = F)
+z.p1 <- (p1-mean(p1))/sd(p1)
+mean(z.p1)
+round(mean(z.p1),10)
+sd(z.p1)
+pnorm(1.8)-pnorm(-1.8)
+hist(z.p1, breaks=100, col=rgb(0,0,0,0))
+abline(v=c(m.p1, -1.8, 1.8), col='red')
+-(pnorm(1.8)-pnorm(-1.8))
+pnorm(1)-pnorm(-1)
 -(pnorm(-1)*2)
-pnorm(-2)
+pnorm(2)-pnorm(-2)
-pnorm(2, lower.tail = F)
 -(pnorm(-2)*2)
-pnorm(-3)
-pnorm(3, lower.tail = F)
 -(pnorm(-3)*2)
+#
+hist(p1, breaks=100, col=rgb(1,1,1,1))
+abline(v=mean(p1),lwd=2)
+abline(v=mean(p1)-sd(p1), lwd=2)
+abline(v=mean(p1)+sd(p1), lwd=2)
+abline(v=c(m.p1-2*sd.p1, m.p1+2*sd.p1), lwd=2, col='red')
+abline(v=c(m.p1-3*sd.p1, m.p1+3*sd.p1), lwd=2, col='green')
 # 68%
@@ Line 36: / Line 116: @@
 c(a, b)
 c(-1, 1)
 # 95%
 c <- qnorm(.05/2)
@@ Line 48: / Line 129: @@
 pnorm(b)-pnorm(a)
+c(a, b)
 pnorm(d)-pnorm(c)
+c(c, d)
 pnorm(f)-pnorm(e)
+c(e, f)
-################################
+qnorm(.5)
-N.p <- 1000000
+qnorm(1)
-m.p <- 100
-sd.p <- 10
-p1 <- rnorm2(N.p, m.p, sd.p)
-mean(p1)
-sd(p1)
-p2 <- rnorm2(N.p, m.p+20, sd.p)
-mean(p2)
-sd(p2)
+################################
 hist(p1, breaks=50, col = rgb(1, 0, 0, 0.5),
      main = "histogram of p1 and p2",)
@@ Line 136: / Line 212: @@
 mean(p1)
 sd(means)
+var(p1)
+sqrt(var(p1)/s.size)
+se.z
+sd(means)
+se.s
 se.z
@@ Line 159: / Line 241: @@
      col = rgb(1, 0, 0, .5))
 abline(v=mean(means), col="black", lwd=2)
-abline(v=mean(p1), col="red", lwd=2)
 # 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"),
        lwd=2)
-abline(v=c(loz1, hiz1, loz2, hiz2, loz3, hiz3),
-       col=c("red","red", "red", "red", "red", "red"),
-       lwd=2)
 round(c(lo1, hi1))
@@ Line 177: / Line 254: @@
 round(c(loz3, hiz3))
-m.sample.i.got <- 101.8
+m.sample.i.got <- mean(means)+ 1.5*sd(means)
+m.sample.i.got
 hist(means,
@@ Line 189: / Line 267: @@
 # m.sample.i.got?
 m.sample.i.got
-pnorm(m.sample.i.got, mean(p1), sd(means), lower.tail = F)
+pnorm(m.sample.i.got, mean(means), sd(means), lower.tail = F)
+pnorm(m.sample.i.got, mean(p1), se.z, lower.tail = F)
 # then, what is the probabilty of getting
@@ Line 199: / Line 278: @@
 abline(v=tmp, col='green', lwd=3)
 * pnorm(m.sample.i.got, mean(p1), sd(means), lower.tail = F)
+m.sample.i.got
 ### one more time
@@ Line 227: / Line 307: @@
 abline(v=tmp, col='green', lwd=2)
 * pnorm(m.sample.i.got, mean(p1), sd(means), lower.tail = F)
 </code>
 ====== output ======
+<WRAP group>
+<WRAP column half>
 <code>
 > rm(list=ls())
 > rnorm2 <- function(n,mean,sd){
 +   mean+sd*scale(rnorm(n))
-+ }
->
-> se <- function(sample) {
-+   sd(sample)/sqrt(length(sample))
 + }
 >
@@ Line 244: / Line 324: @@
 +   sum((x-mean(x))^2)
 + }
-> ################################
->
-> # reminder.
-> pnorm(-1)
-[1] 0.1586553
-> pnorm(1, lower.tail = F)
-[1] 0.1586553
-> 1-(pnorm(-1) + pnorm(1, lower.tail = F))
-[1] 0.6826895
-> 1-(pnorm(-1)*2)
-[1] 0.6826895
->
-> pnorm(-2)
-[1] 0.02275013
-> pnorm(2, lower.tail = F)
-[1] 0.02275013
-> 1-(pnorm(-2)*2)
-[1] 0.9544997
->
-> pnorm(-3)
-[1] 0.001349898
-> pnorm(3, lower.tail = F)
-[1] 0.001349898
-> 1-(pnorm(-3)*2)
-[1] 0.9973002
->
->
-> # 68%
-> a <- qnorm(.32/2)
-> b <- qnorm(1-.32/2)
-> c(a, b)
-[1] -0.9944579  0.9944579
-> c(-1, 1)
-[1] -1  1
-> # 95%
-> c <- qnorm(.05/2)
-> d <- qnorm(1-.05/2)
-> c(c, d)
-[1] -1.959964  1.959964
-> c(-2,2)
-[1] -2  2
-> # 99%
-> e <- qnorm(.01/2)
-> f <- qnorm(1-.01/2)
-> c(e,f)
-[1] -2.575829  2.575829
-> c(-3,3)
-[1] -3  3
->
-> pnorm(b)-pnorm(a)
-[1] 0.68
-> pnorm(d)-pnorm(c)
-[1] 0.95
-> pnorm(f)-pnorm(e)
-[1] 0.99
 >
+</code>
+</WRAP>
+<WRAP column half>
+....
+필요한 펑션
+  * rnorm2는 평균과 표준편차값을 (mean, sd) 갖는 n 개의 샘플을 랜덤하게 추출하는 것
+  * ss는 sum of square 값을 구하는 펑션
+  * ss는 ss/n-1 의 variance 혹은 mean square값을 구할 때의 분자부분으로
+  * error 혹은 residual의 제곱의 합이라고 부를 수 있다 ([[:variance]] 참조)
+</WRAP>
+</WRAP>
+<WRAP group>
+<WRAP column half>
+<code>
 > ################################
 > N.p <- 1000000
@@ Line 317: / Line 358: @@
 [1] 10
 >
-> hist(p1, breaks=50, col = rgb(1, 0, 0, 0.5),
+> m.p1 <- mean(p1)
-+      main = "histogram of p1",)
+> sd.p1 <- sd(p1)
-> abline(v=mean(p1), col="black", lwd=3)
+> var(p1)
+     [,1]
+[1,]  100
 >
+</code>
+</WRAP>
+<WRAP column half>
+  * rnorm2 펑션을 이용하여 p1과 p2 모집단을 구한다.
+  * 각 모집단의 평균은 100과 120이고, 표준편차는 모두 10이다.
+  * 모집단의 크기는 각각 백만이다.
+</WRAP>
+</WRAP>
+===== pnorm =====
+모집단 P1의 히스토그램 (histgram of P1)
+{{:r:pasted:20250910-065721.png}}
+<WRAP group>
+<WRAP column half>
+<code>
+> hist(p1, breaks=100, col=rgb(1,1,1,1))
+> abline(v=mean(p1),lwd=2)
+> abline(v=mean(p1)-sd(p1), lwd=2)
+> abline(v=mean(p1)+sd(p1), lwd=2)
+> abline(v=c(m.p1-2*sd.p1, m.p1+2*sd.p1), lwd=2, col='red')
+> abline(v=c(m.p1-3*sd.p1, m.p1+3*sd.p1), lwd=2, col='green')
 >
-> s.size <- 10
+> # area bet black = 68%
+> # between red = 95%
+> # between green = 99%
 >
-> means.temp <- c()
+> pnorm(m.p1+sd.p1, m.p1, sd.p1)
-> s1 <- sample(p1, s.size, replace = T)
+[1] 0.8413447
-> mean(s1)
+> pnorm(m.p1-sd.p1, m.p1, sd.p1)
-[1] 94.66982
+[1] 0.1586553
-> means.temp <- append(means.temp, mean(s1))
+> pnorm(m.p1+sd.p1, m.p1, sd.p1) -
-> means.temp <- append(means.temp, mean(sample(p1, s.size, replace = T)))
++   pnorm(m.p1-sd.p1, m.p1, sd.p1)
-> means.temp <- append(means.temp, mean(sample(p1, s.size, replace = T)))
+[1] 0.6826895
-> means.temp <- append(means.temp, mean(sample(p1, s.size, replace = T)))
-> means.temp
-[1]  94.66982 100.08316  99.93752  99.83882
 >
-> iter <- 1000000
+> pnorm(m.p1+2*sd.p1, m.p1, sd.p1) -
-> means <- c()
++   pnorm(m.p1-2*sd.p1, m.p1, sd.p1)
-> means <- rep(NA, iter)
+[1] 0.9544997
-> for (i in 1:iter) {
-+   # means <- append(means, mean(sample(p1, s.size, replace = T)))
-+   means[i] <- mean(sample(p1, s.size, replace = T))
-+ }
 >
-> mean(means)
+> pnorm(m.p1+3*sd.p1, m.p1, sd.p1) -
-[1] 100.0042
++   pnorm(m.p1-3*sd.p1, m.p1, sd.p1)
-> sd(means)
+[1] 0.9973002
-[1] 3.163687
-> se.s <- sd(means)
 >
-> lo1 <- mean(means)-se.s
+> m.p1
-> hi1 <- mean(means)+se.s
-> lo2 <- mean(means)-2*se.s
-> hi2 <- mean(means)+2*se.s
-> 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=c(lo1, hi1, lo2, hi2, lo3, hi3),
-+        col=c("green","green", "blue", "blue", "orange", "orange"),
-+        lwd=2)
->
-> # meanwhile . . . .
-> se.s
-[1] 3.163687
->
-> se.z <- sqrt(var(p1)/s.size)
-> se.z <- c(se.z)
-> se.z
-[1] 3.162278
->
-> # 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:
-> # mean(means) = mean(p1)
-> # and
-> # sd(means) = sd(p1) / sqrt(s.size)
-> # that is, sd(means) = se.z
-> # This is called CLT (Central Limit Theorem)
-> mean(means)
-[1] 100.0042
-> mean(p1)
 [1] 100
-> sd(means)
+> sd.p1
-[1] 3.163687
+[1] 10
-> se.z
-[1] 3.162278
 >
-> # because CLT
+> (m.p1-m.p1)/sd.p1
-> loz1 <- mean(p1)-se.z
+[1] 0
-> hiz1 <- mean(p1)+se.z
+> ((m.p1-sd.p1) - m.p1) / sd.p1
-> loz2 <- mean(p1)-2*se.z
+[1] -1
-> hiz2 <- mean(p1)+2*se.z
+> (120-100)/10
-> loz3 <- mean(p1)-3*se.z
+[1] 2
-> hiz3 <- mean(p1)+3*se.z
+> pnorm(1)-pnorm(-1)
+[1] 0.6826895
+> pnorm(2)-pnorm(-2)
+[1] 0.9544997
+> pnorm(3)-pnorm(-3)
+[1] 0.9973002
 >
-> c(lo1, loz1)
+> 1-pnorm(-2)*2
-[1] 96.84055 96.83772
+[1] 0.9544997
-> c(lo2, loz2)
+> pnorm(2)-pnorm(-2)
-[1] 93.67686 93.67544
+[1] 0.9544997
-> c(lo3, loz3)
-[1] 90.51317 90.51317
 >
-> c(hi1, hiz1)
+> pnorm(120, 100, 10)
-[1] 103.1679 103.1623
+[1] 0.9772499
-> c(hi2, hiz2)
+> pnorm(2)-pnorm(-2)
-[1] 106.3316 106.3246
+[1] 0.9544997
-> c(hi3, hiz3)
-[1] 109.4953 109.4868
 >
+>
+>
+>
 >
-> 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)
+</code>
-> abline(v=mean(p1), col="red", lwd=2)
+</WRAP>
-> # abline(v=mean(p2), colo='darkgreen', lwd=3)
+<WRAP column half>
-> abline(v=c(lo1, hi1, lo2, hi2, lo3, hi3),
+....
-+        col=c("green","green", "blue", "blue", "orange", "orange"),
+pnorm 펑션
-+        lwd=2)
+  * 위의 histogram에서 검정색의 선은 p1의 standard deviation값인 10씩 좌우로 그린것 (90과 110)
-> abline(v=c(loz1, hiz1, loz2, hiz2, loz3, hiz3),
+  * 붉은 색선은 80, 120
-+        col=c("red","red", "red", "red", "red", "red"),
+  * 녹색선은 70, 130 선을 말한다
-+        lwd=2)
+  * 고등학교 때, Normal distribution (정상분포의) 경우
+  * 평균을 중심으로 위아래로 SD값 하나씩 간 간격의 probability는 (면적은) 68%
+  * 두 개씩 간 값은 95%
+  * 세 개씩 간 값은 99%  라고 배웠다.
+  * pnorm은 percentage를 구하는 R의 명령어
+  * pnorm(m.p1+sd.p1, m.p1, sd.p1) 은
+    * 정규분포곡선에서
+    * 평균값과 표준편차값을 더한 값 (100+10=110) 을 기준으로
+    * 왼쪽의 부분이 몇 퍼센트인가를 구해주는 명령어
+    * m.p1+sd.p1은 110이므로
+    * 오른 쪽 검정색을 기준으로 왼쪽의 퍼센티지를 묻는 것
+    * 답은 0.8413447
+    * pnorm(m.p1-sd.p1, m.p1, sd.p1) - pnorm(90, 100, 10)이므로
+    * 왼쪽 검정색 선을 기준으로 왼쪽의 퍼센티지를 묻는 것
+    * 답은 0.1586553
+    * 전자에서 후자를 빼 준 값은 두개의 검정색 선의 안쪽 면적으로 묻는 것
+    * 답은 0.6826895
+  * 이 0.6826895 이 우리가 배운 68%
+  * 그렇다면 두 개씩 간 면적은 (붉은 색 선의 안쪽 부분)
+    * pnorm(m.p1+2*sd.p1, m.p1, sd.p1) - pnorm(m.p1-2*sd.p1, m.p1, sd.p1) = 0.9544997
+    * 0.9544997 혹은 95.44997%
+    * 이것이 우리가 배운 95%
+  * 마찬가지로 녹색선 가운데 부분은 pnorm(m.p1+3*sd.p1, m.p1, sd.p1) - pnorm(m.p1-3*sd.p1, m.p1, sd.p1) 로 구할 수 있는데 답은
+  * 0.9973002
+  * pnorm명령어는 pnorm(score, mean, sd)와 같이 사용하여 퍼센티지값을 구할 수 있는데
+  * pnorm(score)만으로 구하면, mean과 sd가 각각 0, 1을 default값으로 하고 생략이 된 것
+  * 위의 p1 모집단도 두번째 방법으로 구해 볼 수 있는데, 그렇게 하기 위해서
+  * p1을 표준점수화 하면 됨
+  * 표준점수화 한다는 뜻은 p1집합의 평균과 표준편차 값을 0과 1로 만든다는 것
+  * p1 모든 원소를 표준점수화하는 것은 각 점수와 평균점수 간의 차이에 sd가 몇개나 들어가는 지 구하는 것. 즉,
+    * $ \dfrac {(\text{score} - \text{m.p1})} {\text{sd.p1}} $
+    * 100점은 표준점수로 0이 된다 $ (100-100) / 10 = 0 $
+    * 110점은 1
+    * 90점은 -1
+    * 115점은 1.5
+    * 130점은 30/10 = 3 점으로 계산될 수 있다.
+  * 그렇다면 위의 pnorm(m.p1+sd.p1, m.p1, sd.p1)은 pnorm(1)과 같은 것
+  * pnorm(110, 100, 10) 이고 이 때 110점은 표준점수로 1 이다.
+  * 따라서 위의 68%, 95%, 99%는 pnorm(1)-pnorm(-1)과 같이 구할 수 있다.
+</WRAP>
+</WRAP>
+===== z score, 표준점수 =====
+<WRAP group>
+<WRAP column half>
+<code>
+> zscore <- (120-100)/10
+> pnorm(zscore)-pnorm(-zscore)
+[1] 0.9544997
+> zscore
+[1] 2
 >
 >
-> round(c(lo1, hi1))
+> pnorm(-1)
-[1]  97 103
+[1] 0.1586553
-> round(c(lo2, hi2))
+> pnorm(1, 0, 1, lower.tail = F)
-[1]  94 106
+[1] 0.1586553
-> round(c(lo3, hi3))
+> pnorm(110, 100, 10, lower.tail = F)
-[1]  91 109
+[1] 0.1586553
+> zscore <- (110-100)/10
+> pnorm(zscore, lower.tail = F)
+[1] 0.1586553
 >
-> round(c(loz1, hiz1))
+> pnorm(118, 100, 10, lower.tail = F)
-[1]  97 103
+[1] 0.03593032
-> round(c(loz2, hiz2))
+> pnorm(18/10, lower.tail = F)
-[1]  94 106
+[1] 0.03593032
-> round(c(loz3, hiz3))
-[1]  91 109
 >
-> m.sample.i.got <- 101.8
+</code>
+</WRAP>
+<WRAP column half>
+....
+</WRAP>
+</WRAP>
+{{:r:pasted:20250910-073339.png}}
+<WRAP group>
+<WRAP column half>
+<code>
+> z.p1 <- (p1-mean(p1))/sd(p1)
+> mean(z.p1)
+[1] -2.319195e-17
+> round(mean(z.p1),10)
+[1] 0
+> sd(z.p1)
+[1] 1
+> pnorm(1.8)-pnorm(-1.8)
+[1] 0.9281394
 >
-> hist(means,
+> hist(z.p1, breaks=100, col=rgb(0,0,0,0))
-+      xlim = c(mean(means)-10*sd(means), mean(means)+10*sd(means)),
+> abline(v=c(m.p1, -1.8, 1.8), col='red')
-+      col = rgb(1, 0, 0, .5))
+> 1-(pnorm(1.8)-pnorm(-1.8))
-> abline(v=mean(means), col="black", lwd=3)
+[1] 0.07186064
-> abline(v=m.sample.i.got, col='darkgreen', lwd=3)
 >
-> # what is the probablity of getting
-> # greater than
-> # m.sample.i.got?
-> m.sample.i.got
-[1] 101.8
-> pnorm(m.sample.i.got, mean(p1), sd(means), lower.tail = F)
-[1] 0.2846929
 >
-> # then, what is the probabilty of getting
+</code>
-> # greater than m.sample.i.got and
+</WRAP>
-> # less than corresponding value, which is
+<WRAP column half>
-> # mean(means) - m.sample.i.got - mean(means)
+....
-> # (green line)
+p1 모집단의 모든 원소를 표준점수화 하기
-> tmp <- mean(means) - (m.sample.i.got - mean(means))
+  * z.p1 <- (p1-mean(p1))/sd(p1)
-> abline(v=tmp, col='green', lwd=3)
+  * 평균과 표준편차는 각각 0, 1이 된다 (mean(z.p1), sd(z.p1))
-> 2 * pnorm(m.sample.i.got, mean(p1), sd(means), lower.tail = F)
+  * 그렇다면 이 표준점수화한 분포에서 -1.8과 1.8 사이를 제외한 바깥 쪽 부분의 면적은 (probability는)
-[1] 0.5693857
+  * 1-(p(1.8)-p(-1.8))과 같이 구할 수 있다.
->
+  * 답은 약 7% 정도
-> ### one more time
+</WRAP>
-> mean(p2)
+</WRAP>
-[1] 120
-> sd(p2)
-[1] 10
-> sample(p2, s.size)
- [1] 112.3761 123.0103 124.5613 116.8093 113.7091 113.2657 126.8168 125.4273 123.8021 118.1890
-> m.sample.i.got <- mean(sample(p2, s.size))
-> m.sample.i.got
-[1] 115.9053
->
-> hist(means,
-+      xlim = c(mean(means)-15*sd(means), mean(means)+15*sd(means)),
-+      col = rgb(1, 0, 0, .5))
-> abline(v=mean(means), col="black", lwd=2)
-> abline(v=m.sample.i.got, col='darkgreen', lwd=2)
->
-> # what is the probablity of getting
-> # greater than
-> # m.sample.i.got?
-> m.sample.i.got
-[1] 115.9053
-> pnorm(m.sample.i.got, mean(p1), sd(means), lower.tail = F)
-[1] 2.484948e-07
->
-> # then, what is the probabilty of getting
-> # greater than m.sample.i.got and
-> # less than corresponding value, which is
-> # mean(means) - m.sample.i.got - mean(means)
-> # (green line)
-> tmp <- mean(means) - (m.sample.i.got - mean(means))
-> abline(v=tmp, col='green', lwd=2)
-> 2 * pnorm(m.sample.i.got, mean(p1), sd(means), lower.tail = F)
-[1] 4.969897e-07
-> 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, breaks=50, col=rgb(0,0,1,.5))
-> p1 <- rnorm2(N.p, m.p, sd.p)
-> mean(p1)
-[1] 100
-> sd(p1)
-[1] 10
-> p2 <- rnorm2(N.p, m.p+20, sd.p)
-> mean(p2)
-[1] 120
-> sd(p2)
-[1] 10
-> 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, breaks=50, col=rgb(0,0,1,.5))
-> 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))
-> s.size <- 10
-> abline(v=mean(p2), col="violet", lwd=3)
-> s.size <- 10
-> means.temp <- c()
-> # means.temp <- c()
-> s1 <- sample(p1, s.size, replace = T)
-> mean(s1)
-[1] 98.98075
-> means.temp <- append(means.temp, mean(s1))
-> means.temp <- c()
-> s1 <- sample(p1, s.size, replace = T)
-> mean(s1)
-[1] 95.23863
-> means.temp <- append(means.temp, mean(s1))
-> means.temp <- append(means.temp, mean(sample(p1, s.size, replace = T)))
-> means.temp <- append(means.temp, mean(sample(p1, s.size, replace = T)))
-> means.temp <- append(means.temp, mean(sample(p1, s.size, replace = T)))
-> means.temp
-[1]  95.23863  98.50240 104.56303  89.11020
-> means.temp <- append(means.temp, mean(sample(p1, s.size, replace = T)))
-> means.temp
-[1]  95.23863  98.50240 104.56303  89.11020  99.51111
-> iter <- 1000000
-> # means <- c()
-> means <- rep(NA, iter)
-> for (i in 1:iter) {
-+   # means <- append(means, mean(sample(p1, s.size, replace = T)))
-+   means[i] <- mean(sample(p1, s.size, replace = T))
-+ }
->
->
-> length(means)
-[1] 1000000
-> length(means)
-[1] 1000000
-> mean(means)
-[1] 99.99929
-> sd(means)
-[1] 3.162556
-> se.s <- sd(means)
-> hist(means, breaks=100, rgb=(.4,.4,.4,.4))
-에러: 예기치 않은 ','입니다 in "hist(means, breaks=100, rgb=(.4,"
-> hist(means, breaks=100, rgb=(.4, .4, .4, .4))
-에러: 예기치 않은 ','입니다 in "hist(means, breaks=100, rgb=(.4,"
-> hist(means, breaks=100, col=rgb(.1, 0, 0, .5))
+<WRAP group>
-> # now we want to get sd of this distribution
+<WRAP column half>
-> lo1 <- mean(means)-se.s
+<code>
-> hi1 <- mean(means)+se.s
+> pnorm(1)-pnorm(-1)
-> lo2 <- mean(means)-2*se.s
-> hi2 <- mean(means)+2*se.s
-> lo3 <- mean(means)-3*se.s
-> hi3 <- mean(means)+3*se.s
-> abline(v=mean(means), col="black", lwd=2)
-> abline(v=mean(means), col="red", lwd=2)
-> 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=c(lo1, hi1, lo2, hi2, lo3, hi3),
-+        col=c("green","green", "blue", "blue", "orange", "orange"),
-+        lwd=2)
-> # meanwhile . . . .
-> se.s
-[1] 3.162556
-> se.z <- sqrt(var(p1)/s.size)
-> se.z <- c(se.z)
-> se.z
-[1] 3.162278
-> # when iter value goes to
-> # unlimited value:
-> # mean(means) = mean(p1)
-> # and
-> # sd(means) = sd(p1) / sqrt(s.size)
-> # that is, sd(means) = se.z
-> # This is called CLT (Central Limit Theorem)
-> mean(means)
-[1] 99.99929
-> rm(list=ls())
-> rnorm2 <- function(n,mean,sd){
-+   mean+sd*scale(rnorm(n))
-+ }
->
-> se <- function(sample) {
-+   sd(sample)/sqrt(length(sample))
-+ }
->
-> ss <- function(x) {
-+   sum((x-mean(x))^2)
-+ }
-> ################################
->
-> # reminder.
-> pnorm(-1)
-[1] 0.1586553
-> pnorm(1, lower.tail = F)
-[1] 0.1586553
-> 1-(pnorm(-1) + pnorm(1, lower.tail = F))
 [1] 0.6826895
 > 1-(pnorm(-1)*2)
 [1] 0.6826895
 >
-> pnorm(-2)
+> pnorm(2)-pnorm(-2)
-[1] 0.02275013
+[1] 0.9544997
-> pnorm(2, lower.tail = F)
-[1] 0.02275013
 > 1-(pnorm(-2)*2)
 [1] 0.9544997
 >
-> pnorm(-3)
-[1] 0.001349898
-> pnorm(3, lower.tail = F)
-[1] 0.001349898
 > 1-(pnorm(-3)*2)
 [1] 0.9973002
 >
+</code>
+</WRAP>
+<WRAP column half>
+....
+</WRAP>
+</WRAP>
+===== qnorm =====
+{{:r:pasted:20250910-075030.png}}
+<WRAP group>
+<WRAP column half>
+<code>
+> #
+> hist(p1, breaks=100, col=rgb(1,1,1,1))
+> abline(v=mean(p1),lwd=2)
+> abline(v=mean(p1)-sd(p1), lwd=2)
+> abline(v=mean(p1)+sd(p1), lwd=2)
+> abline(v=c(m.p1-2*sd.p1, m.p1+2*sd.p1), lwd=2, col='red')
+> abline(v=c(m.p1-3*sd.p1, m.p1+3*sd.p1), lwd=2, col='green')
 >
 > # 68%
@@ Line 654: / Line 606: @@
 > c(-1, 1)
 [1] -1  1
+>
 > # 95%
 > c <- qnorm(.05/2)
@@ Line 671: / Line 624: @@
 > pnorm(b)-pnorm(a)
 [1] 0.68
+> c(a, b)
+[1] -0.9944579  0.9944579
 > pnorm(d)-pnorm(c)
 [1] 0.95
+> c(c, d)
+[1] -1.959964  1.959964
 > pnorm(f)-pnorm(e)
 [1] 0.99
+> c(e, f)
+[1] -2.575829  2.575829
 >
-> ################################
+> qnorm(.5)
-> N.p <- 1000000
+[1] 0
-> m.p <- 100
+> qnorm(1)
-> sd.p <- 10
+[1] Inf
 >
-> p1 <- rnorm2(N.p, m.p, sd.p)
-> mean(p1)
-[1] 100
-> sd(p1)
-[1] 10
 >
-> p2 <- rnorm2(N.p, m.p+20, sd.p)
+</code>
-> mean(p2)
+</WRAP>
-[1] 120
+<WRAP column half>
-> sd(p2)
+....
-[1] 10
+  * qnorm는 pnorm의 반대값을 구하는 명령어
+  * 히스토그램에서 검정 색 부분의 바깥 쪽 부분은 32%이고 왼 쪽의 것은 이것의 반인 16% 이다.
+  * 이 16% 에 해당하는 표준점수가 무엇인가를 묻는 질문이 ''qnorm(.32/2)'' 혹은 ''qnorm(.32/2, 0, 1)''
+  * 원점수일 경우에 대입해서 물어본다면 ''qnorm(.32/2, 100, 10)''과 같이 하면 된다 (혹은 qnorm(.32/2, m.p1, sd.p1) ).
+  * 위의 명령어에 대한 답은 우리가 가장 쉽게 이해한 방법에 의하면
+  * -1 혹은 90 이렇게 될 것이다.
+  * 하지만 위에서 봤던 것처럼 pnorm(1)-pnorm(-1)이 정확히
+  * 0.6826895 인 것처럼, 우리가 이해한 방법은 정확한 값을 구해주지는 않는다.
+  * 마찬가지로 qnorm(0.05/2)와 qnorm(1-(0.05/2))는
+  * 표준편차 값인 1이 왼쪽으로 두개 내려가고 오른 쪽으로 두개 올라간 값의 quotient값을 구하는 것이므로
+  * 즉, 위의 histogram에서 붉은 색 부분의 바깥 쪽 면적에 해당하는 점수를 물어보는 것이므로
+  * -2와 2라고 대답해도 되지만 정확한 답은 아래와 같다 (1.96).
+<code>
+> qnorm(0.05/2)
+[1] -1.959964
+> qnorm(1-(0.05/2))
+[1] 1.959964
 >
+</code>
+  * 이에 해당하는 원점수 (p1 원소에 해당하는) 값은 80과 120이 아니라 아래와 같을 것이다.
+<code>
+> qnorm(0.05/2, 100, 10)
+[1] 80.40036
+> qnorm(1-(0.05/2), 100, 10)
+[1] 119.5996
+>
+</code>
+</WRAP>
+</WRAP>
+<WRAP group>
+<WRAP column half>
+<code>
+> ################################
 > hist(p1, breaks=50, col = rgb(1, 0, 0, 0.5),
 +      main = "histogram of p1 and p2",)
@@ Line 699: / Line 685: @@
 > abline(v=mean(p2), col="violet", lwd=3)
 >
+</code>
+</WRAP>
+<WRAP column half>
+</WRAP>
+</WRAP>
+===== distribution of sample means =====
+<WRAP group>
+<WRAP column half>
+<code>
 > s.size <- 10
 >
@@ Line 704: / Line 700: @@
 > s1 <- sample(p1, s.size, replace = T)
 > mean(s1)
-[1] 107.1742
+[1] 99.64273
 > means.temp <- append(means.temp, mean(s1))
 > means.temp <- append(means.temp, mean(sample(p1, s.size, replace = T)))
@@ Line 711: / Line 707: @@
 > means.temp <- append(means.temp, mean(sample(p1, s.size, replace = T)))
 > means.temp
-[1] 107.17418 101.54022  99.38215  96.19341 100.70208
+[1]  99.64273 107.15516 103.81192 103.12311 105.88372
 >
+</code>
+</WRAP>
+<WRAP column half>
+....
+  * s.size는 10으로 우선 고정하고
+  * 한 샘플을 취하여 그 평균값을 means.temp 메모리에 add한다 (저장하는 것이 아니라 means.temp에 붙힌다).
+  * 그 다음 샘플의 평균을 다시 means.temp에 저장한다
+  * 그 다음 . . . 모두 5번을 하고 출력을 해본다
+</WRAP>
+</WRAP>
+<WRAP group>
+<WRAP column half>
+<code>
 > iter <- 1000000
 > # means <- c()
@@ Line 723: / Line 734: @@
 [1] 1000000
 > mean(means)
-[1] 99.99699
+[1] 99.99544
 > sd(means)
-[1] 3.160973
+[1] 3.161886
 > se.s <- sd(means)
 >
@@ Line 731: / Line 742: @@
 > abline(v=mean(means), col="red", lwd=2)
 >
+</code>
+</WRAP>
+<WRAP column half>
+....
+</WRAP>
+</WRAP>
+{{:r:pasted:20250910-080946.png}}
+<WRAP group>
+<WRAP column half>
+<code>
 > # now we want to get sd of this distribution
 > lo1 <- mean(means)-se.s
@@ Line 748: / Line 771: @@
 +        lwd=2)
 >
+</code>
+</WRAP>
+<WRAP column half>
+....
+  * 위 백만개의 샘플평균이 모인 집합의 히스토그램을 그리고
+  * 그 집합의 표준편차값을 수직선으로 표시하기 위해서
+  * mean(means) +- se.s 와 같은 방법을 쓴 후 그래프로 그린다.
+  * 아래에서 선 하나씩의 길이는 means 집합의 (distribution of sample means)
+  * 표준편차값이다.
+  * 이 표준편차 값을 위에서 sd(means)로 구한 후에 se.s로 저장한 적이 있다.
+  * 그리고 그 값은 3.161886 이었다.
+<code>
+> se.s
+[1] 3.161886
+</code>
+</WRAP>
+</WRAP>
+{{:r:pasted:20250910-081358.png}}
+<WRAP group>
+<WRAP column half>
+<code>
 > # meanwhile . . . .
 > se.s
-[1] 3.160973
+[1] 3.161886
 >
 > se.z <- sqrt(var(p1)/s.size)
@@ Line 758: / Line 805: @@
 >
 > # sd of sample means (sd(means))
-> # is sqrt(var(s1)/s.size) or
-> # sd(s1) / sqrt(s.size)
 > # = se.s
 >
@@ Line 770: / Line 815: @@
 > # This is called CLT (Central Limit Theorem)
 > mean(means)
-[1] 99.99699
+[1] 99.99544
 > mean(p1)
 [1] 100
 > sd(means)
-[1] 3.160973
+[1] 3.161886
+> var(p1)
+     [,1]
+[1,]  100
+> sqrt(var(p1)/s.size)
+         [,1]
+[1,] 3.162278
+> se.z
+[1] 3.162278
+>
+> sd(means)
+[1] 3.161886
+> se.s
+[1] 3.161886
 > se.z
 [1] 3.162278
@@ Line 787: / Line 845: @@
 >
 > c(lo1, loz1)
-[1] 96.83602 96.83772
+[1] 96.83356 96.83772
 > c(lo2, loz2)
-[1] 93.67504 93.67544
+[1] 93.67167 93.67544
 > c(lo3, loz3)
-[1] 90.51407 90.51317
+[1] 90.50978 90.51317
 >
 > c(hi1, hiz1)
-[1] 103.1580 103.1623
+[1] 103.1573 103.1623
 > c(hi2, hiz2)
-[1] 106.3189 106.3246
+[1] 106.3192 106.3246
 > c(hi3, hiz3)
-[1] 109.4799 109.4868
+[1] 109.4811 109.4868
 >
 >
+</code>
+</WRAP>
+<WRAP column half>
+....
+  * 그런데 이 값은 (se.s = 3.161886)
+  * se.z 를 구하는 방법과 거의 같은 값을 갖는다 3.162278
+<code>
+> se.z <- sqrt(var(p1)/s.size)
+> se.z <- c(se.z)
+> se.z
+[1] 3.162278
+</code>
+  * 사실, 우리가 백만 번의 샘플을 취해서 구한 means 집합의 평균과 표준편차 값은
+  * 만약에 백만 번이 아니라 무한 대로 더 큰 숫자를 사용한다고 하면
+  * 위의 se.z 값을 구하는 식의 값을 갖게 된다. 이것을 말로 풀어서 설명하면
+  * 샘플평균들의 집합에서 표준편차 값은
+  * 원래 모집단의 분산값을 샘플사이즈로 나누어준 값에 제곱근을 씌워서 구할 수 있다이다.
+  * 즉, 샘플평균을 모은 집합의 분산값은 그 샘플이 추출된 원래 population의 분산값을 샘플크기로 (sample size) 나누어 준 값이다.
+  * 즉, '' var(means) = var(p1) / s.size ''
+  * 따라서, ''std(means) = sqrt(var(p1) / s.size) ''
+  * 더하여 그 샘플평균 집합의 평균 값은 population의 평균값이 된다
+  * 즉, '' mean(means) = mean(p1) ''
+  * 따라서 lo, hi에 해당하는 means분포의 값을 mean(means) +- sd(means)로 구했었는데,
+  * 샘플평균의 분포를 무한대 번을 했다고 하면 사실 이 값은
+  * mean(p1) +- se.z 로 구하는 것이 정확할 것이다.
+  * 여기서 ''se.z = sqrt(var(p1)/s.size))''
+  * loz1 - hiz1, loz2 - hiz2 값들은 이렇게 구한 값들이다.
+참고
+<code>
+> lo1 <- mean(means)-se.s
+> hi1 <- mean(means)+se.s
+> lo2 <- mean(means)-2*se.s
+> hi2 <- mean(means)+2*se.s
+> lo3 <- mean(means)-3*se.s
+> hi3 <- mean(means)+3*se.s
+</code>
+</WRAP>
+</WRAP>
+{{:r:pasted:20250910-083710.png}}
+<WRAP group>
+<WRAP column half>
+<code>
 > 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(p1), col="red", lwd=2)
 > # 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"),
 +        lwd=2)
-> abline(v=c(loz1, hiz1, loz2, hiz2, loz3, hiz3),
-+        col=c("red","red", "red", "red", "red", "red"),
-+        lwd=2)
->
 >
 > round(c(lo1, hi1))
@@ Line 829: / Line 930: @@
 [1]  91 109
 >
-> m.sample.i.got <- 101.8
+</code>
+</WRAP>
+<WRAP column half>
+....
+</WRAP>
+</WRAP>
+{{:r:pasted:20250910-084125.png}}
+<WRAP group>
+<WRAP column half>
+<code>
+> m.sample.i.got <- mean(means)+ 1.5*sd(means)
+> m.sample.i.got
+[1] 104.7383
 >
 > hist(means,
@@ Line 841: / Line 955: @@
 > # m.sample.i.got?
 > m.sample.i.got
-[1] 101.8
+[1] 104.7383
-> pnorm(m.sample.i.got, mean(p1), sd(means), lower.tail = F)
+> pnorm(m.sample.i.got, mean(means), sd(means), lower.tail = F)
-[1] 0.2845272
+[1] 0.0668072
+> pnorm(m.sample.i.got, mean(p1), se.z, lower.tail = F)
+[1] 0.06701819
 >
 > # then, what is the probabilty of getting
@@ Line 853: / Line 969: @@
 > abline(v=tmp, col='green', lwd=3)
 > 2 * pnorm(m.sample.i.got, mean(p1), sd(means), lower.tail = F)
-[1] 0.5690543
+[1] 0.1339882
+> m.sample.i.got
+[1] 104.7383
 >
+</code>
+</WRAP>
+<WRAP column half>
+....
+  * 만약에 내가 한 샘플을 취해서 평균값을 살펴보니
+  * m.sample.i.got 값이었다고 하자 (104.7383).
+  * 이 값보다 큰 값이거나 아니면
+  * 이 값에 해당하는 평균 반대편 값보다 작은 값이 값이
+  * 나올 확률은 무엇인가?
+  * 즉, 녹색선과 연두색 선 바깥 쪽 부분의 probability 값은?
+  * 아래처럼 구해서 13.4% 정도가 된다
+<code>
+> 2 * pnorm(m.sample.i.got, mean(p1), sd(means), lower.tail = F)
+[1] 0.1339882
+</code>
+</WRAP>
+</WRAP>
+<WRAP group>
+<WRAP column half>
+<code>
 > ### one more time
 > mean(p2)
@@ Line 861: / Line 1001: @@
 [1] 10
 > sample(p2, s.size)
- [1] 115.5752 109.5564 137.8459 115.3346 130.8225 129.1909 118.8858 107.1550 118.8443 109.5731
+ [1] 120.9639 119.8341 134.0344 122.5446 121.4557 113.6820 114.7603 105.2138 122.7027 125.4131
 > m.sample.i.got <- mean(sample(p2, s.size))
 > m.sample.i.got
-[1] 122.4168
+[1] 120.2492
 >
 > hist(means,
@@ Line 876: / Line 1016: @@
 > # m.sample.i.got?
 > m.sample.i.got
-[1] 122.4168
+[1] 120.2492
 > pnorm(m.sample.i.got, mean(p1), sd(means), lower.tail = F)
-[1] 6.621728e-13
+[1] 7.560352e-11
 >
 > # then, what is the probabilty of getting
@@ Line 888: / Line 1028: @@
 > abline(v=tmp, col='green', lwd=2)
 > 2 * pnorm(m.sample.i.got, mean(p1), sd(means), lower.tail = F)
-[1] 1.324346e-12
+[1] 1.51207e-10
->
->
 >
 </code>
+</WRAP>
+<WRAP column half>
+....
+</WRAP>
+</WRAP>