> 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" ) > > > > > > > comment * msrs 의 min(msrs)값을 찾는다 24.975 * 이 최소값이 어느 위치에 있는지 (몇번째 자리에 있는지) 찾는다 which(msrs == min(msrs)) * 이 위치가 151이다 * 이 151번째 사용된 (최소값인 msr값 = mean(x-v)^2)을 결과한) vs값을 찾는다 (vs[151]) * 이 값을 출력하니 50 이고 이 값은 x의 평균값이다. * 이를 plot으로 출력한 것 {{:pasted:20250904-173050.png}}