COMMunication
RESearch.NET

This page is read only. You can view the source, but not change it. Ask your administrator if you think this is wrong.
====== Why n-1 ======
문제.

우리는 모집단의 (population) 평균값을 알고 있다면 샘플의 분산값을 다음과 같이 구할 수 있다. 그리고, 이것을 모집단의 분산값으로 추정할 수 있다. 

\begin{eqnarray*}
\hat{\sigma}^{2} = \dfrac {\displaystyle\sum_{i=1}^{n}{(X_{i}-\mu)}} {n} 
\end{eqnarray*}

그러나, 현재 우리가 가지고 있는 것은 샘플 밖에 없다. 즉, 모집단의 평균은 알지 못하는 상태이기에 모집단 분산을 추정하는 계산에 사용할 수 없다. 따라서 샘플의 평균을 사용한다. 그런데, 샘플의 평균을 사용할 때는 분모에 N 대신에 n-1을 사용해야 한다. 왜 n-1을 사용하는것이 모집단의 분산값 추정에 도움이 되는가가 문제이다. 

\begin{eqnarray*}
\hat{\sigma}^{2} \neq \frac {\displaystyle\sum_{i=1}^{n}{(X_{i}-\overline{X})}} {n} 
\end{eqnarray*}

\begin{eqnarray*}
\hat{\sigma}^{2} = \frac {\displaystyle\sum_{i=1}^{n}{(X_{i}-\overline{X})}} {n-1} 
\end{eqnarray*}

이 모집단의 분산값을 ($\sigma^2$) 대표함을 알아보는 것이 문제이다.


====== 직관적 이해 ======
분산은 $SS/df$ 라고 배웠는데, SS = Sum of Something Square라고 설명하고, 여기서 Something 은 error 값이라고 (개인의 점수를 평균으로 추측했을 때 틀린 만큼의 에러값) 하였다. 그리고 평균으로 예측하는 방법이 가장 작은 오차를 갖는 (예측이 들 틀리는) 방법이라고 하였다 (개인점수 중에서 평균이 제일 많이 나오므로, 평균으로 개인점수를 예측하면 제일 들 틀린다). 따라서 어느 한 집합에서 (샘플에서) 개인점수에서 평균을 빼고 이른 제곱하여 모두 더한 값은 (평균이 아닌) 다른 값으로 구한 값보다 항상 작게 된다 (최소값을 갖는다).

 
이를 그림으로도 설명할 수 있다. 아래에서 녹색의 세로선은 모집단의 평균값이고, 붉은색의 세로선은 3개로 이루어진 샘플의 평균값이다. 그리고 녹색 가로선은 3개의 샘플요소와 모집단평균과의 ($\mu$) 차이값들이고, 적색가로선은 3개의 샘플요소와 샘플평균과의 ($\overline{X}$) 차이값이다. 이 차이값들을 모아서 길이를 비교한 것이 그래프의 하단이다. 적색가로선 세개의 합이 녹색가로선 세개의 합보다 작다. 이는 샘플평균을 사용했을 때와 모집단의 평균을 사용했을 때를 비교하는 것이지만 모집단 평균외에 다른 값을 썼어도 마찬가지이다. 
{{:pasted:20200412-002825.png?800}}

====== 실험적, 수학적 이해 ======
\begin{eqnarray*}
\sum_{i=1}^{n}{(X_{i}-\mu)}  \gt \sum_{i=1}^{n}{(X_{i}-\overline{X})}
\end{eqnarray*}
를 수학적으로 이해하는 방법이다. 우선 실험을 통해서 원하는 것이 무엇인가를 설명한다. 우선 R에서 평균이 20인 (sd = 4) 모집단을 만든다.

<code>
## population parameter 지정
n.p <- 10000
mean.p <- 20
sd.p <- 4
set.seed(23)
p <- rnorm(n.p, mean=mean.p, sd=sd.p)
p <- round(p)
hist(p, freq=F)
curve(dnorm(x, mean=mean(p), sd=sd(p)), add=TRUE, col="blue")
abline(v=mean.p,lwd=3,lty=2, col="red")
</code>
  * 모집단에서 평균이 23인 4개의 원소를 샘플로  취한다. 그리고, 
  * 1부터 40까지의 집합을 만들어 range에 기록해두고 (range = {1,2,3,4,. . .,40}, 이 range에는 위 샘플의 평균인 23이 포함되어 있다) 
  * $\sum{(x-\overline{x})}$ 에서 $\overline(x)$ 대신 1:40 까지의 숫자를 넣어 결과를 구해본다. 즉, SS파트를 구해보는데 샘플의 평균인 23외에 1에서 40까지의 숫자를 대입하여 SS값을 구하여 기록한다는 뜻이다.
  * 이를 plot한다.

<code>
set.seed(1953)
x <-  sample(p, 4)
x
mean(x)

range <- seq(1:40)
ss  <- rep (NA, length(range))
for (i in range) {
    ss[i] <- sum((x-range[i])^2)
}
data <- data.frame(range,ss)
data
min(data$ss) ## ss값이 최소일 때의 x값을 살펴보자 (=mean(x)) = 23
plot(data, lty=1, lwd=1)
abline(v=mean(x),col="red")
</code>

{{:pasted:20200504-202759.png}}
<code>
> n.p <- 10000
> mean.p <- 20
> sd.p <- 4
> set.seed(23)
> p <- rnorm(n.p, mean=mean.p, sd=sd.p)
> p <- round(p)
> head(p)
[1] 21 18 24 27 24 24
> hist(p, freq=F)
> curve(dnorm(x, mean=mean(p), sd=sd(p)), add=TRUE, col="blue")
> abline(v=mean.p,lwd=3,lty=2, col="red")
> 


# 모집단 평균 = 20, sd=4
> set.seed(1953)
> x <-  sample(p, 4)
> x
[1] 27 21 21 23
> mean(x)
[1] 23
> 
> range <- seq(1:40)
> ss  <- rep (NA, length(range))
> for (i in range) {
+     ss[i] <- sum((x-range[i])^2)
+ }
> data <- data.frame(range,ss)
> data
   range   ss
1      1 1960
2      2 1788
3      3 1624
4      4 1468
5      5 1320
6      6 1180
7      7 1048
8      8  924
9      9  808
10    10  700
11    11  600
12    12  508
13    13  424
14    14  348
15    15  280
16    16  220
17    17  168
18    18  124
19    19   88
20    20   60
21    21   40
22    22   28
23    23   24
24    24   28
25    25   40
26    26   60
27    27   88
28    28  124
29    29  168
30    30  220
31    31  280
32    32  348
33    33  424
34    34  508
35    35  600
36    36  700
37    37  808
38    38  924
39    39 1048
40    40 1180
> min(data$ss) ## ss값이 최소일 때의 x값을 살펴보자 (=mean(x)) = 23
[1] 24
> plot(data, lty=1, lwd=1)
> abline(v=mean(x),col="red")

</code>

{{:pasted:20200504-203916.png}}

평균이  20, 표준편차가 4인 집단에서 4개의 샘플을 취하여 그 평균을 구하고, 그 평균을 이용하 SS 부분을 (Sum of Square) 구한다고 했을 때, 평균외에 다른 점수를 이용했을 때 어떻게 되는가를 본 것이다 (range <- seq(1:40)과 같이). ss값이 가장 작았을 때의 x값을 보면 샘플의 평균값임을 알  수 있다.

마지막 그래프에서 가장 작은 기울기값을 갖는 v 값을 구한다고 (derivatives) 가정하고 이해를 하면 수학적으로 이해할 수 있다. ((see https://www.mathsisfun.com/calculus/derivatives-introduction.html))
{{:pasted:20200504-223320.png}}

\begin{eqnarray*}
\dfrac{\text{d}}{\text{dv}} \dfrac{\sum{(x-v)^2}}{n} & = &  \dfrac {\sum{2(x-v)*(-1)}}{n} \\
& = & \dfrac{\sum{-2(x-v)}}{n} \\
& = & -\dfrac{2}{n} \sum{(x-v)} \\
\end{eqnarray*}
위의 식이 0이 (기울기가 0이 되는 부분) 될 때의 v 값을 찾아야 하므로

\begin{eqnarray*}
-\dfrac{2}{n} \sum{(x-v)} & = & 0 \\
\sum{(x-v)} & = & 0 \\
\sum{x} - n*v & = & 0 \\ 
n*v & = & \sum{x} \\
v & = & \dfrac {\sum{x}}{n}  \\
\end{eqnarray*}
위에 따르면, 우리가  찾는 v 값은 샘플의 평균값이 ($\frac {\sum{x}}{n}$) 된다. 따라서, 평균값으로 SS값을 구하게 되면 언제나 가장 작은  값을 취하게 되는 결과를 갖는다. 이렇게 작은 값을 갖는 현상을 보정하려고 n 대신에 n보다 조금 작은 숫자인 n-1을 가지고 SS 부분을 나누어 준다. 

그렇다면  왜 n-2 혹은 n-(1/2)가 아니고  n-1인가? 이를 수학적인 증명을 통해서 살펴보면 다음 장과 같다.





====== 수학적 증명 ======
우선, 

\begin{eqnarray*}
Var[X] & = & E[(X-\mu)^{2}] \\
       & = & E[(X^{2} - 2 X \mu + \mu^{2})] \\
& = & E[X^{2}] - 2 \mu E[X] + E[\mu^2] \\
& = & E[X^{2}] - 2 \mu E[X] + E[\mu^{2}], \;\; \text{because}\; E[X] = \mu \text{, } \; E[\mu^2] = \mu^2, \\
& = & E[X^{2}] - 2 \mu^{2} + \mu^{2}   \\
& = & E[X^{2}] - \mu^{2}
\end{eqnarray*}

이므로

\begin{align}
E\left[X^2\right] & = Var\left[X\right] + \mu^2 \nonumber \\ 
& = \sigma^{2} + \mu^2 \\
\end{align}

마찬가지로 
\begin{align}
Var \left[ \overline{X}\right] & =  E \left[\overline{X}^2 \right] - \left[E(\overline{X})\right]^2 \nonumber \\
& = E\left[\overline{X}^{2}\right] - \mu^{2} \nonumber 
\end{align}

따라서
\begin{align}
E\left[\overline{X}^{2}\right]  & = Var\left[\overline{X}\right] + \mu^2 \nonumber \\ 
& = \frac {\sigma^{2}} {n} + \mu^{2} 
\end{align}

참고로 위에서 $Var\left[\overline{X}\right] = \dfrac {\sigma^{2}} {n} $ 에 해당하는 설명은 [[:mean and variance of the sample mean]] 문서를 볼 것.

----
참고로 Expected value (기대값)와 Variance (분산)의 연산에 과한 법칙으로는 (([[:statistical review]]참조))
<WRAP box 450px>
X,Y are Independent variables.

\begin{align*}
E[aX] & = a E[X] \\
E[X+Y] & = E[X] + E[Y] \\
Var[aX] & = a^{\tiny{2}} Var[X] \\
Var[X+Y] & = Var[X] + Var[Y]  
\end{align*}

</WRAP>
----
우리가 알고자 하는 것은 아래의 식이 population의 parameter인 $\sigma^{2}$ 의 값과 같은가이다.
\begin{align*}
E[s^{2}] & = E \left[\frac{\displaystyle\sum_{i=1}^{n}(X_{i}-\overline{X})^{2}}{n-1} \right] \qquad
\cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot\cdot \;\; (a)  \\
& = \sigma^{2} 
\end{align*}

위의 식에서 일부만을 추출해서 먼저 보자.

\begin{align*}
E \left[\sum{(X_{i}-\overline{X})^{2}} \right] & = E \left[\sum(X_{i}^{2}- 2 X_{i} \overline{X} + \overline{X}^{2})\right] \\
& = E \left[ \sum{X_{i}^2} - \sum{2X_{i} \overline{X}} + \sum {\overline{X^2}}  \right]  \\
& = E \left[ \sum{X_{i}^2} - 2 \overline{X} \sum{X_{i}} + \sum{\overline{X^2}}  \right]  \\
& = E \left[ \sum{X_{i}^2} - 2 \overline{X} \sum{X_{i}} + n \overline{X^2} \right]  \\
& = E \left[ \sum{X_{i}^2} - 2 \overline{X} \cdot (n \overline{X}) + n \overline {X^2} \right] \\
& = E \left[ \sum{X_{i}^2} - n \overline{X}^2 \right] \\
& = \sum {E\left(X_{i}^2\right)} - E\left(n\overline{X}^2\right)  \\
& = \sum {E\left(X_{i}^2\right)} - n E\left(\overline{X}^2\right)  \;\;\; \dots\dots\dots\dots\dots (3)
\end{align*}

한 편, 위의 $(1), (2)$에서 

<WRAP box>
\begin{align*}
E\left[X_{i}^{2}\right] & = \sigma^{2} + \mu^2 \;\;\; \dots\dots\dots\dots\dots (1) \\
E\left[\overline{X}^{2}\right] & = \dfrac {\sigma^{2}} {n} + \mu^{2} \;\;\; \dots\dots\dots\dots\dots (2)
\end{align*}
</WRAP>

위의 $(1), (2)$를 $(3)$에 대입해보면

\begin{align*}
E \left[\sum{(X_{i}-\overline{X})^{2}} \right] & = \sum{E\left(X_{i}^{2}\right)} - n E\left(\overline{X}^{2}\right)  \\
& = \sum{\left(\sigma^{2} + \mu^{2}\right)} - n \left(\dfrac{\sigma^2}{n} + \mu^2\right) \\
& = n\sigma^{2} + n\mu^{2} - \sigma^{2} - n\mu^{2} \\
& = \left(n-1\right) \sigma^{2} 
\end{align*}

위는 식 (a)의 일부이므로 이를 온전한 식에 대입해보면, 
\begin{eqnarray*}
E \left[\sum{(X_{i}-\overline{X})^{2}} \right] & = & (n-1) \sigma^{2} \\
\end{eqnarray*}

\begin{eqnarray*}
E[s^{2}] & = & E \left[ \frac{\displaystyle\sum_{i=1}^{n}(X_{i}-\overline{X})^{2}}{n-1} \right] \\
& = & \dfrac{1}{n-1} E \left[\sum{(X_{i}-\overline{X})^{2}} \right] \\
& = & \dfrac{1}{n-1} (n-1) \sigma^{2} \\
& = & \sigma^{2}
\end{eqnarray*}

그러므로, **n-1로 나눠 준 샘플분산의 (sample's variance) 기대값은**
\begin{eqnarray*}
E(s^2) = \sigma^{2}
\end{eqnarray*}


----
만약에 우리가 population의 variance를 구하듯이 n을 이용한다고 하면, 

\begin{eqnarray*}
E[s^{2}] & = & E \left[ \frac{\displaystyle\sum_{i=1}^{n}(X_{i}-\overline{X})^{2}} {n} \right], \;\;\; \text{note that we use n instead of n-1} \\
& = & \dfrac{1}{n} E \left[\sum{(X_{i}-\overline{X})^{2}} \right] \\
& = & \dfrac{1}{n} (n-1) \sigma^{2} \\
& = & \left(\dfrac{n-1}{n}\right) \sigma^{2} \\
\end{eqnarray*}

즉, 원래 $\sigma^2$ 값보다 조금 작은 값을 갖게 될 것이다 (이를 biased result라고 한다).


====== R 에서 SS값이 최소인 v값을 찾기 ======
<code>
#library(ggplot2)
#library(ggpmisc)

rm(list=ls())
rnorm2 <- function(n,mean,sd){ 
  mean+sd*scale(rnorm(n)) 
}

# set.seed(191)
nx <- 20
mx <- 50
sdx <- mx * 0.15
x <- rnorm2(nx, mx, sdx)

mean(x)
sd(x)
length(x)
hist(x)

x.span <- seq(from = mean(x)-6*sd(x), 
              to = mean(x)+6*sd(x), 
              by = 0.1)

residuals <- function(x, v) {
  return(x - v)
}

ssr <- function(x, v) {
  residuals <- (x - v)
  return(sum(residuals^2))
}

msr <- function(x, v) {
  residuals <- (x - v)
#  return((sum(residuals^2))/(length(x)-1))
  return((mean(residuals^2)))
}

srs <- c() # sum of residuals
ssrs <- c() # sum of square residuals
msrs <- c() # mean square residuals = variance
vs <- c() # the value of v in (x - v)

for (i in x.span) {
  res.x <- residuals(x,i)
  srs.x <- sum(res.x)  
  ssr.x <- ssr(x,i)
  msr.x <- msr(x,i)
  srs <- append(srs, srs.x)
  ssrs <- append(ssrs, ssr.x)
  msrs <- append(msrs, msr.x)
  vs <- append(vs, i)
}
plot(srs)
plot(msrs)
plot(srs)

min(msrs)
min.pos.msrs <- which(msrs == min(msrs))
min.pos.msrs
print(vs[min.pos.msrs])

plot(vs, msrs)
plot(vs, srs)



# the above no gradient
# mse 값으로 계산 rather than sse 
# 후자는 값이 너무 커짐

gradient <- function(x, v){
  residuals = x - v
  dx = -2 * mean(residuals)
  return(list("ds" = dx))
} # function returns ds value

residuals <- function(x, v) {
  return(x - v)
}

ssr <- function(x, v) {
  residuals <- (x - v)
  return(sum(residuals^2))
}

msr <- function(x, v) {
  residuals <- (x - v)
  return((sum(residuals^2))/(length(x)-1))
#  return(mean(residuals^2))
}

# pick one random v in (x-v)
v <- rnorm(1)
# Train the model with scaled features
learning.rate = 1e-1

ssrs <- c()
msrs <- c()
mres <- c()
vs <- c()
# Record Loss for each epoch:
zx <- (x-mean(x))/sd(x)

nlen <- 100
for (epoch in 1:nlen) {
  residual <- residuals(zx, v)
  ssr.x <- ssr(zx, v)
  msr.x <- msr(zx, v)
  ssrs <- append(ssrs, ssr.x)
  msrs <- append(msrs, msr.x)
  
  grad <- gradient(zx, v)
  
  step.v <- grad$ds * learning.rate
  v <- v - step.v
  vs <- append(vs, v)
}

tail(srs)
tail(msrs)
tail(ssrs)
tail(vs)

plot(srs)
plot(msrs)
plot(ssrs)
plot(vs)
# scaled
v
v.orig <- (v*sd(x))+mean(x) 
v.orig

</code>
====== output ======
<code>
> #library(ggplot2)
> #library(ggpmisc)
> 
> rm(list=ls())
> rnorm2 <- function(n,mean,sd){ 
+   mean+sd*scale(rnorm(n)) 
+ }
> 
> # set.seed(191)
> nx <- 20
> mx <- 50
> sdx <- mx * 0.15
> x <- rnorm2(nx, mx, sdx)
> 
> mean(x)
[1] 50
> sd(x)
[1] 7.5
> length(x)
[1] 20
> hist(x)
> 
> x.span <- seq(from = mean(x)-6*sd(x), 
+               to = mean(x)+6*sd(x), 
+               by = 0.1)
> 
> residuals <- function(x, v) {
+   return(x - v)
+ }
> 
> ssr <- function(x, v) {
+   residuals <- (x - v)
+   return(sum(residuals^2))
+ }
> 
> msr <- function(x, v) {
+   residuals <- (x - v)
+ #  return((sum(residuals^2))/(length(x)-1))
+   return((mean(residuals^2)))
+ }
> 
> srs <- c() # sum of residuals
> ssrs <- c() # sum of square residuals
> msrs <- c() # mean square residuals = variance
> vs <- c() # the value of v in (x - v)
> 
> for (i in x.span) {
+   res.x <- residuals(x,i)
+   srs.x <- sum(res.x)  
+   ssr.x <- ssr(x,i)
+   msr.x <- msr(x,i)
+   srs <- append(srs, srs.x)
+   ssrs <- append(ssrs, ssr.x)
+   msrs <- append(msrs, msr.x)
+   vs <- append(vs, i)
+ }
> plot(srs)
> plot(msrs)
> plot(srs)
> 
> min(msrs)
[1] 53.4375
> min.pos.msrs <- which(msrs == min(msrs))
> min.pos.msrs
[1] 451
> print(vs[min.pos.msrs])
[1] 50
> 
> plot(vs, msrs)
> plot(vs, srs)
> 
> 
> 
> # the above no gradient
> # mse 값으로 계산 rather than sse 
> # 후자는 값이 너무 커짐
> 
> gradient <- function(x, v){
+   residuals = x - v
+   dx = -2 * mean(residuals)
+   return(list("ds" = dx))
+ } # function returns ds value
> 
> residuals <- function(x, v) {
+   return(x - v)
+ }
> 
> ssr <- function(x, v) {
+   residuals <- (x - v)
+   return(sum(residuals^2))
+ }
> 
> msr <- function(x, v) {
+   residuals <- (x - v)
+   return((sum(residuals^2))/(length(x)-1))
+ #  return(mean(residuals^2))
+ }
> 
> # pick one random v in (x-v)
> v <- rnorm(1)
> # Train the model with scaled features
> learning.rate = 1e-1
> 
> ssrs <- c()
> msrs <- c()
> mres <- c()
> vs <- c()
> # Record Loss for each epoch:
> zx <- (x-mean(x))/sd(x)
> 
> nlen <- 100
> for (epoch in 1:nlen) {
+   residual <- residuals(zx, v)
+   ssr.x <- ssr(zx, v)
+   msr.x <- msr(zx, v)
+   ssrs <- append(ssrs, ssr.x)
+   msrs <- append(msrs, msr.x)
+   
+   grad <- gradient(zx, v)
+   
+   step.v <- grad$ds * learning.rate
+   v <- v - step.v
+   vs <- append(vs, v)
+ }
> 
> tail(srs)
[1] -890 -892 -894 -896 -898 -900
> tail(msrs)
[1] 1 1 1 1 1 1
> tail(ssrs)
[1] 19 19 19 19 19 19
> tail(vs)
[1] 2.936258e-11 2.349006e-11 1.879204e-11 1.503363e-11 1.202690e-11 9.621523e-12
> 
> plot(srs)
> plot(msrs)
> plot(ssrs)
> plot(vs)
> # scaled
> v
[1] 9.621523e-12
> v.orig <- (v*sd(x))+mean(x) 
> v.orig
[1] 50
> 
> 
</code>
====== tags ======
{{tag>"research methods" "조사방법론" "표준편차" "n-1" "자유도" "degrees of freedom" "n-1" "표준오차"}}