Table of Contents
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$) 대표함을 알아보는 것이 문제이다.
직관적 이해
분산은 $ \text{SS}/ \text{df} $ 라고 배웠다. SS = Sum of Something Square라고 설명하고, 여기서 Something 은 error (개인의 점수를 평균으로 추측했을 때 틀린 만큼의 에러값) 하였다. 이렇게 평균을 가지고 개인점수를 예측하는 것이 가장 작은 오차를 갖는 방법이라고 하였다 (개인점수 중에서 평균이 제일 많이 나오므로, 평균으로 개인점수를 예측하면 제일 들 틀린다). 따라서 어느 한 집합에서 (샘플에서) 개인점수에서 평균을 빼고 이를 제곱하여 모두 더한 값은 (평균이 아닌) 다른 값으로 구한 값보다 항상 작게 된다 (최소값을 갖는다).
이를 그림으로도 설명할 수 있다. 아래에서 녹색의 세로선은 모집단의 평균값이고, 붉은색의 세로선은 3개로 이루어진 샘플의 평균값이다. 그리고 녹색 가로선은 3개의 샘플요소와 모집단평균과의 ($\mu$) 차이값들이고, 적색가로선은 3개의 샘플요소와 샘플평균과의 ($\overline{X}$) 차이값이다. 이 차이값들을 모아서 길이를 비교한 것이 그래프의 하단이다. 적색가로선 세개의 합이 녹색가로선 세개의 합보다 작다. 이는 샘플평균을 사용했을 때와 모집단의 평균을 사용했을 때를 비교하는 것이지만 모집단 평균외에 다른 값을 썼어도 마찬가지이다.
실험적, R에서 시뮬레이션으로 이해
아래 output의 코멘트를 읽을 것
rm(list=ls()) rnorm2 <- function(n,mean,sd){ mean+sd*scale(rnorm(n)) } # set.seed(191) nx <- 1000 mx <- 50 sdx <- mx * 0.1 sdx # 5 x <- rnorm2(nx, mx, sdx) # x <- rnorm2(1000, 50, 5) 와 동일 mean(x) sd(x) length(x) hist(x) x.span <- seq(from = mean(x)-3*sd(x), to = mean(x)+3*sd(x), by = .1) x.span residuals <- function(x, v) { return(x - v) } # sum of square residual 값을 # 구하는 펑션 ssr <- function(x, v) { residuals <- (x - v) return(sum(residuals^2)) } # mean square residual 값을 # 구하는 펑션 (mean square # residual = variance) msr <- function(x, v) { residuals <- (x - v) return((mean(residuals^2))) } ssrs <- c() # sum of square residuals msrs <- c() # mean square residuals = variance vs <- c() # the value of v in (x - v) # x.span의 값들을 v값으로 삼아 sum(x-x.span)^2 처럼 구하면 # SS값을 구한 것이 된다. 우리가 배운 SS값은 x.span의 값으로 # 샘플의 평균을 사용했을 때의 residual 값이다. x.span은 # 샘플의 평균을 중심으로 여러가지 값을 사용하는 것을 가정한다. for (i in x.span) { res.x <- residuals(x,i) msr.x <- msr(x,i) msrs <- append(msrs, msr.x) vs <- append(vs, i) } # 아래 plot은 SS값들이나 (두번째는) MS값들을 v값이 변화함에 # 따라서 (x.span의 범위에 따라서 변화) 어떻게 변화하는지 # 구한 것 plot(msrs) # v값이 x.span에 따라서 변화하여 대입되었을 때의 # MS값들을 (msr 펑션으로 구한 mean square값) # 모아 놓은 값이 msrs msrs # 아래는 위에서 계산한 msr 값들을 저장한 msrs값들 중에서 최소값이 # 되는 것을 찾은 것. 우리는 이것이 샘플의 평균임을 안다. min(msrs) # 최소값일 때의 위치 (msrs에서 몇번째인지) min.pos.msrs <- which(msrs == min(msrs)) min.pos.msrs # msr 최소값이 구해졌을 때 사용된 v값 vs[min.pos.msrs]
output
> rm(list=ls()) > > rnorm2 <- function(n,mean,sd){ + mean+sd*scale(rnorm(n)) + } > > # set.seed(191) > nx <- 1000 > mx <- 50 > sdx <- mx * 0.1 > sdx # 5 [1] 5 > x <- rnorm2(nx, mx, sdx) > # x <- rnorm2(1000, 50, 5) 와 동일 > > mean(x) [1] 50 > sd(x) [1] 5 > length(x) [1] 1000 > hist(x) >
SS = sum(x-mean(x))^2 인데, mean(x)을 즉, x집합의 평균을, x 원소값을 예측하는데 (빼는데) 사용하면 SS값이 최소값이 된다고 하였다. 이것을 R에서 simulation으로 알아보기 위해서 mean(x) 대신에 다른 숫자들을 넣어보려고 한다. 이를 v라고 하면 sum(x-v)^2이라는 SS값들을 구해서 비교하려는 것이다. 대입할 숫자들은 (v) mean(x) +- 3 sd(x) 를 범위로 하고, 그 범위의 시작 숫자에서 (시작은 mean(x)-3sd(x)가 된다) 0.1씩 증가시키면서 대입하고, 각 숫자마다 (처음 숫자는 35, 다음 숫자는 35.1 . . . ) SS값을 구해서 저장하여 그것을 그래프로 그려보고 최소값이 어떤 것인지 보는 것이 진행하려는 작업이다.
단, 이 코드에서 SS대신 MS값을 (SS값을 n으로 나눈 값, 즉, variance값 혹은 Mean Square값) 구해서 보려고 하는데 이것은 같은 의미를 갖는다. 즉, 모든 SS값들에 n을 공토으로 나누어준 값을 저장하고 비교하려는 것이다.
> x.span <- seq(from = mean(x)-3*sd(x), + to = mean(x)+3*sd(x), + by = .1) > x.span [1] 35.0 35.1 35.2 35.3 35.4 35.5 35.6 35.7 35.8 [10] 35.9 36.0 36.1 36.2 36.3 36.4 36.5 36.6 36.7 [19] 36.8 36.9 37.0 37.1 37.2 37.3 37.4 37.5 37.6 [28] 37.7 37.8 37.9 38.0 38.1 38.2 38.3 38.4 38.5 [37] 38.6 38.7 38.8 38.9 39.0 39.1 39.2 39.3 39.4 [46] 39.5 39.6 39.7 39.8 39.9 40.0 40.1 40.2 40.3 [55] 40.4 40.5 40.6 40.7 40.8 40.9 41.0 41.1 41.2 [64] 41.3 41.4 41.5 41.6 41.7 41.8 41.9 42.0 42.1 [73] 42.2 42.3 42.4 42.5 42.6 42.7 42.8 42.9 43.0 [82] 43.1 43.2 43.3 43.4 43.5 43.6 43.7 43.8 43.9 [91] 44.0 44.1 44.2 44.3 44.4 44.5 44.6 44.7 44.8 [100] 44.9 45.0 45.1 45.2 45.3 45.4 45.5 45.6 45.7 [109] 45.8 45.9 46.0 46.1 46.2 46.3 46.4 46.5 46.6 [118] 46.7 46.8 46.9 47.0 47.1 47.2 47.3 47.4 47.5 [127] 47.6 47.7 47.8 47.9 48.0 48.1 48.2 48.3 48.4 [136] 48.5 48.6 48.7 48.8 48.9 49.0 49.1 49.2 49.3 [145] 49.4 49.5 49.6 49.7 49.8 49.9 50.0 50.1 50.2 [154] 50.3 50.4 50.5 50.6 50.7 50.8 50.9 51.0 51.1 [163] 51.2 51.3 51.4 51.5 51.6 51.7 51.8 51.9 52.0 [172] 52.1 52.2 52.3 52.4 52.5 52.6 52.7 52.8 52.9 [181] 53.0 53.1 53.2 53.3 53.4 53.5 53.6 53.7 53.8 [190] 53.9 54.0 54.1 54.2 54.3 54.4 54.5 54.6 54.7 [199] 54.8 54.9 55.0 55.1 55.2 55.3 55.4 55.5 55.6 [208] 55.7 55.8 55.9 56.0 56.1 56.2 56.3 56.4 56.5 [217] 56.6 56.7 56.8 56.9 57.0 57.1 57.2 57.3 57.4 [226] 57.5 57.6 57.7 57.8 57.9 58.0 58.1 58.2 58.3 [235] 58.4 58.5 58.6 58.7 58.8 58.9 59.0 59.1 59.2 [244] 59.3 59.4 59.5 59.6 59.7 59.8 59.9 60.0 60.1 [253] 60.2 60.3 60.4 60.5 60.6 60.7 60.8 60.9 61.0 [262] 61.1 61.2 61.3 61.4 61.5 61.6 61.7 61.8 61.9 [271] 62.0 62.1 62.2 62.3 62.4 62.5 62.6 62.7 62.8 [280] 62.9 63.0 63.1 63.2 63.3 63.4 63.5 63.6 63.7 [289] 63.8 63.9 64.0 64.1 64.2 64.3 64.4 64.5 64.6 [298] 64.7 64.8 64.9 65.0
x-mean(x) = residual = error
sum(residual^2) = SS (sum of square)
SS/n = variance, mean square (ms, MS)
이 residual을 구하기 위해서 쓰는 mean(x)의 대체값들을 (v값들) x.span에 모아 놓은 것이다.
이 값을 출력해보았는데 35.1 에서 시작하여 65에서 끝나며, 0.1씩 증가한다.
> > residuals <- function(x, v) { + return(x - v) + } > > # sum of square residual 값을 > # 구하는 펑션 > ssr <- function(x, v) { + residuals <- (x - v) + return(sum(residuals^2)) + } > > # mean square residual 값을 > # 구하는 펑션 (mean square > # residual = variance) > msr <- function(x, v) { + residuals <- (x - v) + return((mean(residuals^2))) + } >
- 이 후 쓸 function들. (x-v) = residual이라고 부르니까 이 residual을 모으는 function
- function ssr = x 집합과 v값 (x.span의 한 숫자)를 인수를 주었을 때 구할 수 있는 Sum of Square값들 (실제로는 사용하지 않는다. 대신 msr 펑션으로 MS값을 구한다).
- function msr = 의 ssr을 n으로 나누어 구한 mean square residual을 (분산) 구하는 function
> ssrs <- c() # sum of square residuals > msrs <- c() # mean square residuals = variance > vs <- c() # the value of v in (x - v) > > # x.span의 값들을 v값으로 삼아 sum(x-x.span)^2 처럼 구하면 > # SS값을 구한 것이 된다. 우리가 배운 SS값은 x.span의 값으로 > # 샘플의 평균을 사용했을 때의 residual 값이다. x.span은 > # 샘플의 평균을 중심으로 여러가지 값을 사용하는 것을 가정한다. > > for (i in x.span) { + res.x <- residuals(x,i) + msr.x <- msr(x,i) + msrs <- append(msrs, msr.x) + vs <- append(vs, i) + } > # 아래 plot은 SS값들이나 (두번째는) MS값들을 v값이 변화함에 > # 따라서 (x.span의 범위에 따라서 변화) 어떻게 변화하는지 > # 구한 것 > > plot(vs, msrs) >
comment
- x.span의 처음값인 35.1을 넣어서 (x-v)를 구한 후
- msr 펑션으로 mean square residual 값을 구한다.
- 그리고 이 값을 어딘가에 (msrs) 저장한 후
- 그 다음 value인 35.2값을 v 대신에 넣어 다시
- msr값을 구하여 위의 msrs에 추가하여 저장한다.
- 이것을 x.span의 모든값에 걸쳐 진행하면
- msrs에는 x.span의 모든 값을 대응하여 구한
- msr값들이 저장된다. 이 msr값 중에서 최소값을 찾고
- 이 최소값을 구할 때 쓴 v값을 (x.span 중 하나의 값) 찾고자 한다.
- 이를 위해서 for 를 이용한 loop문을 쓴다.
> # v값이 x.span에 따라서 변화하여 대입되었을 때의 > # MS값들을 (msr 펑션으로 구한 mean square값) > # 모아 놓은 값이 msrs > msrs [1] 249.975 246.985 244.015 241.065 238.135 [6] 235.225 232.335 229.465 226.615 223.785 [11] 220.975 218.185 215.415 212.665 209.935 [16] 207.225 204.535 201.865 199.215 196.585 [21] 193.975 191.385 188.815 186.265 183.735 [26] 181.225 178.735 176.265 173.815 171.385 [31] 168.975 166.585 164.215 161.865 159.535 [36] 157.225 154.935 152.665 150.415 148.185 [41] 145.975 143.785 141.615 139.465 137.335 [46] 135.225 133.135 131.065 129.015 126.985 [51] 124.975 122.985 121.015 119.065 117.135 [56] 115.225 113.335 111.465 109.615 107.785 [61] 105.975 104.185 102.415 100.665 98.935 [66] 97.225 95.535 93.865 92.215 90.585 [71] 88.975 87.385 85.815 84.265 82.735 [76] 81.225 79.735 78.265 76.815 75.385 [81] 73.975 72.585 71.215 69.865 68.535 [86] 67.225 65.935 64.665 63.415 62.185 [91] 60.975 59.785 58.615 57.465 56.335 [96] 55.225 54.135 53.065 52.015 50.985 [101] 49.975 48.985 48.015 47.065 46.135 [106] 45.225 44.335 43.465 42.615 41.785 [111] 40.975 40.185 39.415 38.665 37.935 [116] 37.225 36.535 35.865 35.215 34.585 [121] 33.975 33.385 32.815 32.265 31.735 [126] 31.225 30.735 30.265 29.815 29.385 [131] 28.975 28.585 28.215 27.865 27.535 [136] 27.225 26.935 26.665 26.415 26.185 [141] 25.975 25.785 25.615 25.465 25.335 [146] 25.225 25.135 25.065 25.015 24.985 [151] 24.975 24.985 25.015 25.065 25.135 [156] 25.225 25.335 25.465 25.615 25.785 [161] 25.975 26.185 26.415 26.665 26.935 [166] 27.225 27.535 27.865 28.215 28.585 [171] 28.975 29.385 29.815 30.265 30.735 [176] 31.225 31.735 32.265 32.815 33.385 [181] 33.975 34.585 35.215 35.865 36.535 [186] 37.225 37.935 38.665 39.415 40.185 [191] 40.975 41.785 42.615 43.465 44.335 [196] 45.225 46.135 47.065 48.015 48.985 [201] 49.975 50.985 52.015 53.065 54.135 [206] 55.225 56.335 57.465 58.615 59.785 [211] 60.975 62.185 63.415 64.665 65.935 [216] 67.225 68.535 69.865 71.215 72.585 [221] 73.975 75.385 76.815 78.265 79.735 [226] 81.225 82.735 84.265 85.815 87.385 [231] 88.975 90.585 92.215 93.865 95.535 [236] 97.225 98.935 100.665 102.415 104.185 [241] 105.975 107.785 109.615 111.465 113.335 [246] 115.225 117.135 119.065 121.015 122.985 [251] 124.975 126.985 129.015 131.065 133.135 [256] 135.225 137.335 139.465 141.615 143.785 [261] 145.975 148.185 150.415 152.665 154.935 [266] 157.225 159.535 161.865 164.215 166.585 [271] 168.975 171.385 173.815 176.265 178.735 [276] 181.225 183.735 186.265 188.815 191.385 [281] 193.975 196.585 199.215 201.865 204.535 [286] 207.225 209.935 212.665 215.415 218.185 [291] 220.975 223.785 226.615 229.465 232.335 [296] 235.225 238.135 241.065 244.015 246.985 [301] 249.975 >
comment
- msrs값에 저장된 msr값들 (mean square residual값들) 중에서
- 가장 작은 값을 찾아서 그 값을 구하도록 한 v값을 찾고자 한다.
- msrs값을 눈으로 살펴보기에는 너무 힘드므로 . . . .
> # 아래는 위에서 계산한 msr 값들을 저장한 msrs값들 중에서 최소값이 되는 것을 찾은 > # 것. 우리는 이것이 샘플의 평균임을 안다. > min(msrs) [1] 24.975 > # 최소값일 때의 위치 (msrs에서 몇번째인지) > min.pos.msrs <- which(msrs == min(msrs)) > min.pos.msrs [1] 151 > # msr 최소값이 구해졌을 때 사용된 v값 > vs[min.pos.msrs] [1] 50 > > plot(vs, msrs, cex=1, lwd=1, lty=3) > abline(v=vs[min.pos.msrs]) > text(x=50, y=150, "msr gets minimal value, when v = 50" ) > > > > > > >
다음으로 MS값을 구하는 식인
- y = sum( (x-v)^2 ) / n 값을
v에 대해서v를 가지고 y를 미분하면 (dy/dv)- 변화하는 v값마다의 기울기 값을 구할 수 있는데
- 이 값이 0이 되는 지점의 v값이 무엇인지를 구하면 위의 R코드에서 구한 값을 구할 수 있게 된다. 아래는 그 과정이다.
- 그래프에서 가장 작은 기울기값을 갖는 v 값을 구한다고 (derivatives) 가정하고 이해를 하면 수학적으로 이해할 수 있다. 1)
\begin{eqnarray*}
\dfrac{dy}{dv} = \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*}
즉, 기울기가 0이 될때의 v값은 x집합의 평균값일 때이다.
R에서 msr이 최소값이 되는 v값 효율적으로 찾기
위의 미분을 이용한 방법을 써서 MS값이 최소값이 되는 순간을 R에서 구현해 볼 수 있다. 아래 아웃풋의 코멘트 부분을 읽을 것.
why n-1 gradient explanation
# the above no gradient gradient <- function(x, v){ residuals = x - v # y = (sum(x-v)^2)/n 혹은 (mean(x-v)^2) # 의 식이 ms값을 구하는 식인데 # 이를 v에 대해서 미분하면 chain rule을 써야 한다 # 자세한 것은 http://commres.net/wiki/estimated_standard_deviation # 문서 중에 미분 부분 참조 # dy/dv = ( 2(x-v)*-1 ) / n chain rule # dy/dv = -2 (x-v) / n = -2 (mean(residual)) dx = -2 * mean(residuals) # return(list("ds" = dx)) return(dx) } # function returns ds value zx <- (x-mean(x))/sd(x) # pick one random v in (x-v) v <- rnorm(1) # Train the model with scaled features learning.rate = 1e-1 msrs <- c() vs <- c() nlen <- 75 for (epoch in 1:nlen) { residual <- residuals(zx, v) msr.x <- msr(zx, v) msrs <- append(msrs, msr.x) grad <- gradient(zx, v) step.v <- grad * learning.rate # v <- v - step.v # 그 다음 v값 vs <- append(vs, v) # v값 저장 } tail(msrs) tail(vs) plot(vs, msrs, type='b') # scaled vs # 변화하는 v 값들의 집합 vs.orig <- (vs*sd(x))+mean(x) vs.orig # 마지막 v값이 최소값에 근접한 값 v v.orig <- (v*sd(x))+mean(x) v.orig plot(vs.orig, msrs, type='b')
output
> # the above no gradient > > gradient <- function(x, v){ + residuals = x - v + # y = (sum(x-v)^2)/n 혹은 (mean(x-v)^2) + # 의 식이 ms값을 구하는 식인데 + # 이를 v에 대해서 미분하면 chain rule을 써야 한다 + # 자세한 것은 http://commres.net/wiki/estimated_standard_deviation + # 문서 중에 미분 부분 참조 + # dy/dv = ( 2(x-v)*-1 ) / n chain rule + # dy/dv = -2 (x-v) / n = -2 (mean(residual)) + dx = -2 * mean(residuals) + # return(list("ds" = dx)) + return(dx) + } # function returns ds value > >
comment
이 R script의 목적은 v값이 최소값이 되는 지점을 자동적으로 찾아보려는 것이다. 이것을 위해서 우선 v값으로 사용할 첫 점수를 랜덤하게 구한 후 (아래 그래프에서 빨간색 지점), 자동적으로 그 다음 v 점수를 찾고 (녹색지점), 그 다음 v 점수를 찾고 (황금색 지점), . . . 이런 과정을 계속하면서 각 v 점수에서의 msr값을 구해서 이에 해당하는 v값을 찾아 보려고 한다. 빨간색, 녹색, 황금색, . . . 이를 자동적으로 구하기 위해서 두가지 방법을 사용하는데 그 것이
- gradient function과
- learning_rate 값이다.
gradient 펑션은 dy/dv 의 연쇄 미분식인 (chain rules) -2(x-v) / n = -2 mean(res) 값을 구하는 것이다. 이렇게 구한 값에 learning_rate값을 곱한후, 이것을 먼저 사용한 v값에서 (빨간색 지점) 빼 주어 다음 v값으로 (녹색지점) 사용하려고 한다. 이 녹색지점에서의 v값을 사용했을 때의 gradient값을 구한 후 다시 이값에 learning_rate인 0.1을 곱하여 그다음 스텝의 값을 얻고, 이 값을 바로 전의 v값에서 빼 준 값을 그 다음 v값으로 사용한다. 이렇게 구하는 v값들은 0.1씩 곱해주는 효과때문에 오른 쪽으로 옮겨가는 지점이 “점진적으로 줄어들게 되고” 이 지점이 msr의 최소값이 되는 지점으로 가게 된다.
> zx <- (x-mean(x))/sd(x) > # pick one random v in (x-v) > v <- rnorm(1) > >
comment
- 랜덤하게 v값을 찾음
v ← rnorm(1)
- 원래는 mean(x)값 근처의 값을 랜덤하게 골라야 하므로
v ← rnorm(1, mean(x), sd(x))
와 같이 써야 하지만 (현재의 경우,rnorm(1, 50, 5)
가 된다)- 그렇게 하질 않고 x 집합의 원소들을 표준점수화 한 후
zx ← (x-mean(x))/sd(x)
(이렇게 표준점수화 하면 x 변인의 평균이 0, 표준편차가 1 이 되는 집합으로 변한다) v ← rnorm(1, 0, 1)
로 구한다. 뒤의 인자인 0, 1 은 default이모로 생략.- 이렇게 하는 이유는 혹시 나중에 다른 x집합에 똑같은 작업을 하더라도 그 집합의 평균과 표준편차를 사용하지 않고
- 단순히
rnorm(1)
을 이용해서 찾으려고 하는 것이다.
> # Train the model with scaled features > learning.rate = 1e-1 >
comment
- 이 0.1은 gradient function으로 구한 해당 v값에 대한 y 미분값을 (기울기 값을) 구한 후, 여기에 곱하기 위해서 지정한다.
> msrs <- c() > vs <- c() > > nlen <- 75 > for (epoch in 1:nlen) { + residual <- residuals(zx, v) + msr.x <- msr(zx, v) + msrs <- append(msrs, msr.x) + + grad <- gradient(zx, v) + step.v <- grad * learning.rate # + v <- v - step.v # 그 다음 v값 + vs <- append(vs, v) # v값 저장 + } > > tail(msrs) [1] 0.999 0.999 0.999 0.999 0.999 0.999 > tail(vs) [1] 6.415945e-08 5.132756e-08 4.106205e-08 3.284964e-08 2.627971e-08 2.102377e-08 > > plot(vs, msrs) > >
comment
- nlen는 그냥 자의적으로 지정한다. 여기서는 75로 했다.
- for 문에서 처음 v값은 위에서 랜덤으로 구한 값이다 (v).
- 이 v값으로 gradient 펑션의 아웃풋 값을 구하고 (-2(residual)) =
grad ← gradient(x, v)
- 이 값에 learning_rate값을 곱한 값을 구하여
step.v ← grad * learaning_rate
- 이 값을 원래 v값에서 빼준 후에
- 다시 (for 문에서 반복하는 동안) v값으로 바꾼 후
v ← v - step.v
- for문을 nlen번 만큼 반복한다.
- 이 과정에서 저장한
- msr값들과
- vs값들의 마지막 6개를 살펴본다.
- v 값으로 x집합의 평균값을 사용하는 것이 최소 msr값이 된다는 것이 맞다면
- v 값은 0이 될것이다 (왜냐하면 x집합을 zx로 바꿨기 때문에, 즉 평균이 0이고 sd값이 1일 집합으로 바꿨기 때문에)
- 아래 그래프의 각 포인트는 v값의 이동을 나타내는데 grad*learning_rate의 영향으로 점진적으로 하가하여 최소값으로 도달한다.
> # scaled > vs # 변화하는 v 값들의 집합 [1] 3.119260e-01 2.495408e-01 1.996326e-01 1.597061e-01 1.277649e-01 1.022119e-01 8.176952e-02 [8] 6.541561e-02 5.233249e-02 4.186599e-02 3.349279e-02 2.679424e-02 2.143539e-02 1.714831e-02 [15] 1.371865e-02 1.097492e-02 8.779935e-03 7.023948e-03 5.619158e-03 4.495327e-03 3.596261e-03 [22] 2.877009e-03 2.301607e-03 1.841286e-03 1.473029e-03 1.178423e-03 9.427383e-04 7.541907e-04 [29] 6.033525e-04 4.826820e-04 3.861456e-04 3.089165e-04 2.471332e-04 1.977066e-04 1.581652e-04 [36] 1.265322e-04 1.012258e-04 8.098061e-05 6.478449e-05 5.182759e-05 4.146207e-05 3.316966e-05 [43] 2.653573e-05 2.122858e-05 1.698286e-05 1.358629e-05 1.086903e-05 8.695226e-06 6.956181e-06 [50] 5.564945e-06 4.451956e-06 3.561565e-06 2.849252e-06 2.279401e-06 1.823521e-06 1.458817e-06 [57] 1.167054e-06 9.336428e-07 7.469143e-07 5.975314e-07 4.780251e-07 3.824201e-07 3.059361e-07 [64] 2.447489e-07 1.957991e-07 1.566393e-07 1.253114e-07 1.002491e-07 8.019931e-08 6.415945e-08 [71] 5.132756e-08 4.106205e-08 3.284964e-08 2.627971e-08 2.102377e-08 > vs.orig <- (vs*sd(x))+mean(x) > vs.orig [1] 51.55963 51.24770 50.99816 50.79853 50.63882 50.51106 50.40885 50.32708 50.26166 50.20933 [11] 50.16746 50.13397 50.10718 50.08574 50.06859 50.05487 50.04390 50.03512 50.02810 50.02248 [21] 50.01798 50.01439 50.01151 50.00921 50.00737 50.00589 50.00471 50.00377 50.00302 50.00241 [31] 50.00193 50.00154 50.00124 50.00099 50.00079 50.00063 50.00051 50.00040 50.00032 50.00026 [41] 50.00021 50.00017 50.00013 50.00011 50.00008 50.00007 50.00005 50.00004 50.00003 50.00003 [51] 50.00002 50.00002 50.00001 50.00001 50.00001 50.00001 50.00001 50.00000 50.00000 50.00000 [61] 50.00000 50.00000 50.00000 50.00000 50.00000 50.00000 50.00000 50.00000 50.00000 50.00000 [71] 50.00000 50.00000 50.00000 50.00000 50.00000 > > # 마지막 v값이 최소값에 근접한 값 > v [1] 2.102377e-08 > v.orig <- (v*sd(x))+mean(x) > v.orig [1] 50 > > plot(vs.orig, msrs, type='b') > > >
그렇다면 왜 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 (분산)의 연산에 과한 법칙으로는 2)
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*}
우리가 알고자 하는 것은 아래의 식이 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)$에서
\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*}
위의 $(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라고 한다).