> 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으로 출력한 것