Cross Validation 및 Bias Variance Trade off 시뮬레이션

Cross Validation

 이번 포스팅에서는 Cross Validation에 대해 다루도록 하겠습니다. Cross Validation이란 Train, Test set을 나누는 과정을 피자조각처럼 나누어 할당하는 것으로 이해하시면 편합니다.

위 그림을 보시면 Data를 4조각으로 나눈 다음에 서로 돌아가면서 Test Set의 역할을 돌아가면서 하는 것을 확인할 수가 있습니다. 즉 쉽게 생각하면 Train set으로 모형을 만들고 Test set으로 검증하는 단계를 쪼갠 만큼 반복하는 것으로 생각하시면 됩니다. 이렇게 하는 이유는 모형의 타당성(Validation)확보에 있습니다. ’이렇게 나누고, 저렇게 나누어도 비슷한 결과가 나온다. 그러니 이 모형은 타당하다.’의 주장을 입증하기 위함입니다. 여기서 마지막으로 계산되는 값은 Cross-Validateion Error(CVE)라고 합니다.

\[ CVE = \frac{1}{n}\sum_{k=1}^{K}n_{k}MSE_{k}\\ MSE_{k}=\frac{1}{n_{k}}\sum_{i\in C_{k}}(y_i-\hat{y}_{i}^{[-k]})^2 \]

library(ggplot2)
library(dplyr)
library(reshape)
library(factoextra)
library(FactoMineR)

AMES = read.csv("D:\\Dropbox\\DATA SET(Dropbox)\\Ames_City.csv",stringsAsFactors = FALSE)

MS.Zoning	Overall.Qual	Year.Built	Year.Remod.Add	Exter.Qual
RL	6	1960	1960	3
RH	5	1961	1961	3
RL	6	1958	1958	3
RL	7	1968	1968	4
RL	5	1997	1998	3
RL	6	1998	1998	3

Exter.Cond	Garage.Qual	Garage.Cond	Sale.Type	Sale.Condition	SalePrice
3	3	3	WD	Normal	215000
3	3	3	WD	Normal	105000
3	3	3	WD	Normal	172000
3	3	3	WD	Normal	244000
3	3	3	WD	Normal	189900
3	3	3	WD	Normal	195500

다음의 데이터로 Cross Validation Error을 계산해보도록 하겠습니다. 조각은 6조각으로 나누도록 하겠습니다. 이때 4-CV가 됩니다.

K = 4

Sample = sample(rep(1:K,length = nrow(AMES)))

MSE_V = c()

for(i in 1:K){
  
  TRAIN = AMES[which(Sample != i),]
  TEST = AMES[which(Sample == i),]
  
  NonLinear = lm(SalePrice ~ poly(Overall.Qual,2), data = TRAIN)
  NonLinear_Predicted = predict(NonLinear, newdata = TEST)
  NonLinear_E = TEST$SalePrice - NonLinear_Predicted
  MSE_V[i] = sum(NonLinear_E^2)/nrow(TEST)
  
}

계산은 매우 간단합니다. 데이터를 4조각으로 나눈 다음에 for문을 활용하여 차례대로 Modeling 및 검증을 반복하면 됩니다.

names(MSE_V) = 1:K
MSE_V

         1          2          3          4 
1872540857 1798261932 2149470897 2055542227 

위 결과값을 확인하면 4번 모두 비슷하면서도 다른 값이 나온 것을 확인할 수가 있습니다. 이제 이 값들에 대해 평균을 취해주면 CVE를 계산할 수가 있습니다. 사실 SalePrice의 값이 매우 크기 때문에 좀 비교하기가 애매하기는 합니다.(사실 이런 문제점들을 개선하기 위하여 데이터에 대해 log 변환, Normalization, Scaling 등을 진행해줍니다. 저도 비교하기가 어렵기 때문에 log변환을 하고 진행하도록 하겠습니다.

K = 4

Sample = sample(rep(1:K,length = nrow(AMES)))

MSE_V = c()

for(i in 1:K){
  
  TRAIN = AMES[which(Sample != i),]
  TEST = AMES[which(Sample == i),]
  
  NonLinear = lm(log(SalePrice) ~ poly(Overall.Qual,2), data = TRAIN)
  NonLinear_Predicted = predict(NonLinear, newdata = TEST)
  NonLinear_E = log(TEST$SalePrice) - NonLinear_Predicted
  MSE_V[i] = sum(NonLinear_E^2)/nrow(TEST)
  
}

names(MSE_V) = 1:K
MSE_V

         1          2          3          4 
0.04370626 0.05303571 0.06134567 0.05373700 

log변환을 하고 회귀분석을 하니 결과가 비슷하게 나오는 것을 확인할 수가 있습니다.

mean(MSE_V)

[1] 0.05295616

CVE는 0.053정도가 나오는 것을 확인할 수가 있습니다.

심화문제

사실 위 문제는 어찌보면 매우 간단합니다. 그러니 이번에는 좀 더 난해한 산점도를 난수 생성한 다음에 Cross Validation을 해보도록 하겠습니다.

set.seed(10000)
X = runif(n = 1000, min = -10, max = 10)
Y = sin(X) + rnorm(n = 1000, mean = 0, sd = 0.3)

ggplot(NULL) +
  geom_point(aes(x = X, y = Y), col = 'royalblue', alpha = 0.5) +
  geom_smooth(aes(x = X, y = Y), se = FALSE, col = 'red') +
  theme_bw() +
  theme(text = element_text(size = 15, face = "bold"))

보기만해도 아찔합니다. 이걸 회귀선으로 한번 적합해보는 연습을 해보도록하겠습니다. 먼저, 이 회귀선에 대해 우리는 몇차항을 적용시켜야될지도 잘 모르겠습니다. 그러니 항차를 몇차항으로 주어야할지도 실험을 통해 확인해야됩니다. 이 문제를 해결하기 위한 for문의 parameter는 다음과 같습니다.

• Replication : R(실험 반복 횟수)
  
• Degree : D(비선형 회귀의 항차)
  
• Kfold : K(K-CV)

어려울 수 있지만 어차피 해야되는 것은 동일합니다. 핵심 알맹이는 어차피 회귀분석을 진행하고 MSE를 계산하는 것입니다. 위에서 만들었던 코드를 기준으로 윗단에는 Replication, 아랫단에는 Degree에 해당되는 for문을 추가시켜주면 됩니다. 다만, for문에서 산출되는 값들을 저장해야될 데이터 공간들은 헷갈릴 수가 있습니다. 최대한 주석을 달도록 하겠습니다.

set.seed(12123)

R = 10 
D = 20
K = 5
d = 5

MSE_R = list()
VAR_R = list()

for(r in 1:R){
  
  Sample = sample(rep(1:K,length = length(X)))
  
  # Cross Validatiton 값 저장공간
  
  MSE_CV = list()
  VAR_CV = list()
  
  for(k in 1:K){
    
    # Degree에 따른 값 저장공간
    
    MSE_Degree = c()
    VAR_Degree = c()
    
    MSE_Degree[1] = NA
    VAR_Degree[1] = NA
    
    TRAIN_X = X[which(Sample != k)]
    TRAIN_Y = Y[which(Sample != k)]
  
    TRAIN_DF = data.frame(
      Y = TRAIN_Y,
      X = TRAIN_X
    )
    
    TEST_X = X[which(Sample == k)]
    TEST_Y = Y[which(Sample == k)]
    
    TEST_DF = data.frame(
      Y = TEST_Y,
      X = TEST_X
    )
    
    for(d in 2:D){
      
      NonLinear = lm(Y ~ poly(X,d), data = TRAIN_DF)
      NonLinear_Predicted = predict(NonLinear, newdata = TEST_DF)
      NonLinear_E = TEST_Y - NonLinear_Predicted
      
      # Degree에 따른 MSE, VAR값 저장
      
      MSE_Degree[d] = sum(NonLinear_E^2)/nrow(TEST)
      VAR_Degree[d] = var(NonLinear_Predicted)
      
    }
    
    # CV에 따른 MSE, VAR값 저장
    
    MSE_CV[[k]] = MSE_Degree
    VAR_CV[[k]] = VAR_Degree
    
  }
  
  
  ### CV DATAFRAME
  
  CVE_DF = data.frame(
    
    DEGREE = 1:D
    
  )
  
  for(i in 1:k){
    
    CVE_DF = cbind(CVE_DF,MSE_CV[[i]])
    
  }
  
  
  ### VAR DATAFRAME
  
  VAR_DF = data.frame(
    
    DEGREE = 1:D
    
  )
  
  for(i in 1:k){
    
    VAR_DF = cbind(VAR_DF,VAR_CV[[i]])
    
  }
  
  ### CVE 계산
  
  CVE = rowMeans(CVE_DF[,2:(k+1)],na.rm = TRUE)
  VAR_E = rowMeans(VAR_DF[,2:(k+1)],na.rm = TRUE)
  
  MSE_R[[r]] = CVE
  VAR_R[[r]] = VAR_E
  
  
}


### 최종저장

CVE_R = data.frame(
  
  DEGREE = 1:D
  
)

for(i in 1:R){
  
  CVE_R = cbind(CVE_R,MSE_R[[i]])
  colnames(CVE_R)[i+1] = paste0("R",i)
  
}



### VAR DATAFRAME

VAR_RDF = data.frame(
  
  DEGREE = 1:D
  
)


for(i in 1:R){
  
  VAR_RDF = cbind(VAR_RDF,VAR_R[[i]])
  colnames(VAR_RDF)[i+1] = paste0("R",i)
  
}

• Replication 및 Degree에 따른 MSE값

DEGREE	R1	R2	R3	R4	R5	R6	R7	R8	R9	R10
1	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN
2	0.1505846	0.1499027	0.1504386	0.1497744	0.1500148	0.1501945	0.1507492	0.1503504	0.1503153	0.1502130
3	0.1482794	0.1468028	0.1472652	0.1468860	0.1469249	0.1471360	0.1476292	0.1474276	0.1473920	0.1470653
4	0.1485622	0.1471725	0.1474680	0.1473158	0.1471559	0.1470297	0.1476818	0.1475707	0.1473677	0.1471298
5	0.1391625	0.1385581	0.1396196	0.1386186	0.1381183	0.1386690	0.1389559	0.1382828	0.1384111	0.1385914
6	0.1390609	0.1388286	0.1396835	0.1387651	0.1380262	0.1386734	0.1389471	0.1383186	0.1383974	0.1387185
7	0.0838003	0.0839175	0.0847815	0.0842244	0.0835188	0.0848734	0.0846038	0.0832664	0.0839124	0.0839613
8	0.0839440	0.0839219	0.0851612	0.0845975	0.0837037	0.0849635	0.0847760	0.0838529	0.0840191	0.0839645
9	0.0319144	0.0320818	0.0319655	0.0323834	0.0318788	0.0319648	0.0320355	0.0318056	0.0317372	0.0318909
10	0.0318910	0.0320625	0.0319033	0.0325826	0.0320069	0.0319252	0.0321572	0.0318513	0.0316846	0.0319690
11	0.0241303	0.0242865	0.0240654	0.0242384	0.0242570	0.0241958	0.0241161	0.0240841	0.0240662	0.0242489
12	0.0242371	0.0242653	0.0240841	0.0243060	0.0242558	0.0242242	0.0241579	0.0240373	0.0240262	0.0242011
13	0.0237422	0.0238904	0.0236658	0.0238945	0.0238118	0.0237294	0.0237270	0.0235902	0.0235129	0.0237641
14	0.0237533	0.0239487	0.0237016	0.0239053	0.0239407	0.0237231	0.0237072	0.0236278	0.0235417	0.0237762
15	0.0237692	0.0239394	0.0236940	0.0238452	0.0239322	0.0237322	0.0236804	0.0236176	0.0235335	0.0237465
16	0.0238169	0.0239327	0.0236669	0.0238212	0.0239393	0.0237082	0.0236565	0.0236243	0.0235526	0.0236812
17	0.0238991	0.0240117	0.0236796	0.0239166	0.0239529	0.0237503	0.0238217	0.0236644	0.0236333	0.0237794
18	0.0240197	0.0240267	0.0237464	0.0239823	0.0239378	0.0237638	0.0238308	0.0237373	0.0236994	0.0237866
19	0.0240571	0.0240714	0.0238329	0.0240346	0.0239750	0.0238298	0.0239189	0.0237851	0.0237767	0.0238071
20	0.0240642	0.0241640	0.0238613	0.0240520	0.0240946	0.0238787	0.0240159	0.0238160	0.0238014	0.0239010

• Replication 및 Degree에 따른 VAR값

DEGREE	R1	R2	R3	R4	R5	R6	R7	R8	R9	R10
1	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN
2	0.0126601	0.0122431	0.0122784	0.0124058	0.0124037	0.0123512	0.0126738	0.0125032	0.0126739	0.0124374
3	0.0260583	0.0243532	0.0247395	0.0252123	0.0251393	0.0250928	0.0248301	0.0253012	0.0254878	0.0248240
4	0.0282069	0.0261047	0.0262737	0.0273439	0.0262403	0.0265485	0.0262172	0.0265422	0.0271956	0.0263800
5	0.0590864	0.0585828	0.0586823	0.0587936	0.0583646	0.0589342	0.0586565	0.0585167	0.0591890	0.0590144
6	0.0597935	0.0593994	0.0592672	0.0598257	0.0593254	0.0594934	0.0594045	0.0593940	0.0596818	0.0600686
7	0.2621596	0.2613105	0.2648010	0.2632597	0.2614296	0.2647483	0.2626203	0.2604441	0.2615481	0.2623707
8	0.2632256	0.2612035	0.2660622	0.2659911	0.2627353	0.2649340	0.2637129	0.2622435	0.2616926	0.2625974
9	0.4424007	0.4448092	0.4460271	0.4486640	0.4469597	0.4437126	0.4458941	0.4459205	0.4450708	0.4460666
10	0.4428475	0.4450725	0.4462037	0.4496585	0.4471969	0.4442821	0.4471311	0.4464641	0.4456042	0.4466714
11	0.4721709	0.4728556	0.4728458	0.4759385	0.4747073	0.4734135	0.4743794	0.4740634	0.4729080	0.4729056
12	0.4727379	0.4728358	0.4728545	0.4770919	0.4741948	0.4737821	0.4736914	0.4742519	0.4730865	0.4731545
13	0.4743411	0.4744297	0.4750917	0.4785824	0.4759438	0.4757865	0.4745462	0.4759252	0.4751324	0.4750693
14	0.4743730	0.4742259	0.4752973	0.4784469	0.4763013	0.4756659	0.4745727	0.4759987	0.4751593	0.4747536
15	0.4749089	0.4747823	0.4754328	0.4788516	0.4766413	0.4758734	0.4751588	0.4761521	0.4751737	0.4750620
16	0.4751531	0.4748010	0.4757108	0.4787704	0.4769315	0.4758477	0.4755685	0.4766293	0.4757013	0.4751752
17	0.4747507	0.4755001	0.4758447	0.4786808	0.4768966	0.4761246	0.4742820	0.4764192	0.4757524	0.4755558
18	0.4747280	0.4757274	0.4763591	0.4785141	0.4768052	0.4760244	0.4741015	0.4765627	0.4764474	0.4755906
19	0.4746639	0.4754801	0.4763730	0.4784033	0.4768960	0.4759772	0.4738536	0.4763712	0.4766403	0.4756408
20	0.4746923	0.4754637	0.4763796	0.4790945	0.4764408	0.4762275	0.4734080	0.4764884	0.4767895	0.4755656

반복횟수, Degree, K-CV를 모두 고려한 실험코드를 완성했습니다. 이제 결과를 확인해보도록 하겠습니다.

CVE_R %>%
  melt(id.vars = c("DEGREE")) %>%
  na.omit() %>%
  ggplot() +
  geom_point(aes(x = DEGREE, y = value, col = variable)) +
  geom_line(aes(x = DEGREE, y = value, col = variable)) +
  labs(col = "Replications") + ylab("MSE(CV)") +
  theme_bw() +
  theme(text = element_text(size = 15, face = "bold"),
        legend.position = "top")

MSE의 경우 10차항 이후부터는 크게 변동이 없는 것을 확인할 수가 있습니다.

VAR_RDF %>%
  melt(id.vars = c("DEGREE")) %>%
  na.omit() %>%
  ggplot() +
  geom_point(aes(x = DEGREE, y = value, col = variable)) +
  geom_line(aes(x = DEGREE, y = value, col = variable)) +
  labs(col = "Replications") + ylab("VAR(CV)") +
  theme_bw() +
  theme(text = element_text(size = 15, face = "bold"),
        legend.position = "top")

Bias-Variance trade-off의 관계에 따라서 Var도 매우 폭등한 것을 확인할 수가 있습니다. 이제 두 그래프를 동시에 비교해보도록 하겠습니다.

rbind(
  CVE_R %>%
    melt(id.vars = c("DEGREE")) %>%
    mutate(G = "MSE") %>%
    na.omit(),
  VAR_RDF %>%
    melt(id.vars = c("DEGREE")) %>%
    mutate(G = "VAR") %>%
    na.omit()
) %>%
  group_by(DEGREE,G) %>%
  summarise(value = mean(value)) %>%
  ggplot() +
  geom_point(aes(x = DEGREE, y = value, col = G),size = 5) +
  geom_line(aes(x = DEGREE, y = value, col = G),size = 1.2) +
  labs(col = "Value") + ylab("MSE & VAR") +
  theme_bw() +
  theme(text = element_text(size = 15, face = "bold"),
        legend.position = "top")

아름다운 Bias-Variance Trade off가 만들어졌습니다. 만약에 CV의 결과만을 확인해보고 싶으신 경우는 MSE_CV, VAR_CV의 데이터를 확인하시면 됩니다.