혼합모형 EM simulation

 

 

 

 

1.

 

시뮬레이션이라는 말은 가끔 들어 봤죠. 그러나 실제로 어떻게 작업하는지 잘 감이 안 올겁니다. 제가 두 집단 혼합모형의 경우 EM 알고리즘이 잘 들어 맞는지 한번 시뮬레이션 해봤습니다. 세 가지를 했는데요. 첫 번째는 초기값을 모수 값에 대한 별 정보가 없는 상태에서 했습니다. 그런데 결과가 아주 안 좋았습니다.

 

두 번째는 생성된 데이터를 밀도추정(density estimation)이라는 것을 해서 초기값을 처음 보다 적절하게 정했습니다. 이 경우 매우 잘 맞췄습니다.

 

세 번째는 첫 번째와 같은데 분포가 식별이 잘 되는 즉, 두 집단이 평균이 첫 번째 보다 차이가 많게끔 했습니다. 그러니까 이 경우도 EM일 잘 맞췄습니다.

 

 

시뮬레이션 상황은 다음과 같습니다. 우리가 관찰한 키 Y의 데이터는 남성 아동과 남성 성인의 두 집단에서 왔다고 가정했습니다. 그리고 아동의 비중은 0.1, 성인의 비중은 0.9로 했습니다. 그리고 아동의 키는 정규분포 N(150, 2), 성인의 키는 정규분포 N(170, 5)로 했습니다. 즉 우리가 하는 일은 오로지 키의 데이터만을 가지고 이 모수 0.1, 0.9, 150, 2, 170, 5 값을 맞춰야 합니다. 데이터는 위의 정보 하에서 100개만 추출했습니다.

 

S, R 코드만 마지막에 제시했습니다. 이 코드는 Davison 책 내용 그대로 따라고요, Davison 책을 첨부합니다. 이 책의 missing value 부분을 보면 EM 부분이 있습니다.

 

davison모델선택책.pdf

 

 

 

2. 첫 번째 시도

 

첫 번째는 초기값을 적당히 넣었습니다. 특별한 정보가 없다고 가정하고 아동의 비중을 0.5, 성인의 비중을 0.5, 그리고 아동과 성인의 평균과 표준편차의 초기값을 넣기 위해 키 데이터를 크기별로 진열하고 반으로 짤라 키가 작은 집단, 키가 큰 집단으로 나눠 여기서 구한 평균과 표준편차를 아동의 키 평균, 표준편차, 성인의 키 평균, 표준편차의 초기값으로 하였습니다. EM 반복 횟수는 50으로 했는데 표에서는 20번까지의 결과만 제시했습니다.

 

	아동 비중	성인 비중	아동 평균	성인 평균	아동 표준편차	성인 표준편차
진짜 값	0.1	0.9	150	170	2	5
EM 반복수
1(초기 값)	0.50	0.50	161.88	173.38	6.95	2.90
2	0.50	0.50	162.72	172.63	7.92	3.08
3	0.50	0.50	163.16	172.14	8.45	3.12
4	0.50	0.50	163.36	171.83	8.75	3.15
5	0.49	0.51	163.44	171.63	8.92	3.19
6	0.48	0.52	163.44	171.51	9.03	3.23
7	0.47	0.53	163.40	171.42	9.11	3.28
8	0.46	0.54	163.34	171.35	9.16	3.33
9	0.46	0.54	163.26	171.31	9.20	3.37
10	0.45	0.55	163.17	171.27	9.23	3.41
11	0.44	0.56	163.08	171.23	9.26	3.45
12	0.43	0.57	162.98	171.20	9.28	3.48
13	0.43	0.57	162.89	171.18	9.30	3.51
14	0.42	0.58	162.80	171.16	9.31	3.54
15	0.42	0.58	162.71	171.13	9.33	3.57
16	0.41	0.59	162.63	171.12	9.34	3.60
17	0.41	0.59	162.55	171.10	9.35	3.62
18	0.40	0.60	162.47	171.08	9.36	3.65
19	0.40	0.60	162.39	171.07	9.36	3.67
20	0.39	0.61	162.31	171.06	9.37	3.69

 

20번 반복했는데도 실제 모수 값을 전혀 못 찾아 가고 있죠.

 

 

 

 

3. 초기값을 적절하게 한 경우

 

진짜 모수 값을 의한 분포 가정에서 데이터 Y를 100개 생성하였죠, 그래서 이걸 밀도추정이라는 것을 해 봤습니다. 그래서 아동 비중을 0.2, 성인 비중을 0.8로 하고 아동의 평균은 작은 키부터 20명의 평균, 그리고 표준편차, 성인은 21명부터 평균과 표준편차를 구해 초기값으로 하였습니다. 아래 밀도추정에서 아동과 성인의 평균값을 어느 정도 알 수 있지만 현실 데이터에서는 이렇게 잘 나오는 경우가 없습니다.

 

 

 

	아동 비중	성인 비중	아동 평균	성인 평균	아동 표준편차	성인 표준편차
진짜 값	0.1	0.9	150	170	2	5
EM 반복수
1(초기 값)	0.20	0.80	155.84	171.06	6.25	4.06
2	0.18	0.82	155.26	170.73	6.52	4.31
3	0.16	0.84	154.54	170.61	6.26	4.39
4	0.15	0.85	153.80	170.54	5.72	4.43
5	0.14	0.86	153.04	170.47	5.02	4.46
6	0.13	0.87	152.33	170.40	4.25	4.50
7	0.12	0.88	151.73	170.33	3.54	4.56
8	0.12	0.88	151.29	170.26	3.00	4.63
9	0.11	0.89	150.95	170.19	2.60	4.70
10	0.11	0.89	150.65	170.12	2.22	4.78
11	0.10	0.90	150.35	170.05	1.77	4.87
12	0.10	0.90	150.16	170.00	1.40	4.92
13	0.10	0.90	150.14	169.99	1.35	4.93
14	0.10	0.90	150.14	169.99	1.35	4.93
15	0.10	0.90	150.14	169.99	1.35	4.93
16	0.10	0.90	150.14	169.99	1.35	4.93
17	0.10	0.90	150.14	169.99	1.35	4.93
18	0.10	0.90	150.14	169.99	1.35	4.93
19	0.10	0.90	150.14	169.99	1.35	4.93
20	0.10	0.90	150.14	169.99	1.35	4.93

 

결과를 보면 10번 정도만 하면 추정값들이 모수와 거의 비슷합니다. 단지 아동의 표준편차만 조금 안 좋게 나옵니다.

 

 

 

 

4. 식별이 잘 되는 경우

 

초기값은 처음처럼 안 좋게 하되 분포 자체가 처음보다 좋게 해서 EM을 적용해 봤습니다. 처음의 경우 아동의 키 평균은 150, 성인의 키 평균이 170으로 가정했는데 이번에는 아동의 키 평균은 150, 성인의 키 평균은 190으로 하였습니다. 이 가정에서 키 데이터 100개를 생성하고 밀도 추정을 하였습니다. 아래 그림입니다.

 

  

  

이 경우 EM 알고리즘 적용 결과입니다. 아동과 성인의 비중의 초기값 0.5, 0.5로 하였습니다.

	아동 비중	성인 비중	아동 평균	성인 평균	아동 표준편차	성인 표준편차
진짜 값	0.1	0.9	150	190	2	5
EM 반복수
1(초기 값)	0.50	0.50	178.01	194.56	14.42	3.91
2	0.49	0.51	179.53	192.88	15.89	4.10
3	0.44	0.56	179.26	191.82	17.06	4.38
4	0.37	0.63	178.13	191.08	18.37	4.66
5	0.30	0.70	176.24	190.62	19.73	4.91
6	0.24	0.76	173.64	190.38	20.85	5.10
7	0.20	0.80	170.70	190.28	21.42	5.23
8	0.18	0.82	167.92	190.25	21.41	5.33
9	0.16	0.84	165.55	190.26	21.01	5.40
10	0.15	0.85	163.57	190.27	20.35	5.46
11	0.14	0.86	161.76	190.29	19.46	5.52
12	0.13	0.87	159.85	190.32	18.19	5.59
13	0.12	0.88	157.50	190.36	16.12	5.67
14	0.11	0.89	154.36	190.39	12.40	5.75
15	0.10	0.90	151.00	190.35	6.32	5.79
16	0.10	0.90	149.89	190.32	1.93	5.81
17	0.10	0.90	149.88	190.32	1.93	5.81
18	0.10	0.90	149.88	190.32	1.93	5.81
19	0.10	0.90	149.88	190.32	1.93	5.81
20	0.10	0.90	149.88	190.32	1.93	5.81

 

결과를 보면 원래 식별이 잘 되는 경우 초기값이 안 좋아도 20회 정도면 충분히 모수값과 비슷하게 추정하다는 것을 알 수 있습니다.

 

 

현실에서는 식별이 잘 안되는 경우도 많이 있을 수 있으니까 사전에 초기값을 조심해서 넣어야 할 겁니다. 변수가 하나면 밀도추정 분포를 그린다든지 변수가 여러개면 군집분석을 실시하면 상당한 좋은 초기값을 구할 수 있을 겁니다.

 

 

 

5. S, R 코드

 

다음은 혼합모형의 S, R 코드입니다. 함수형태로 하지 않았습니다. 함수형태로 하면 범용적으로 사용할 수 있겠죠. 아래 내용을 노트패드에 적어 mixture로 저장하고 S나 R의 명령문 창에서 source("mixture.txt")를 치면 위 결과들이 구해집니다.

truea1<-0.1;

truea2<-0.9;

 

trueu1<-150;

trueu2<-190;

truesigma1<-2;

truesigma2<-5;

n<-100;

y1<-rnorm(0.1*100, trueu1, truesigma1);# data generate

y2<-rnorm(0.9*100, trueu2, truesigma2);

y<-sort(c(y1,y2)); #initial value를 위해 y를 낮은값부터 배열

y1<-y[1:50]; # 아동에게 키가 작은 50명 배당

y2<-y[51:100]; # 성인에게 키가 큰 50명 배당

ema1<-0.5;#a의 초기값으로 각각 0.5

ema2<-0.5;

emu1<-mean(y1);

emu2<-mean(y2);

emsigma1<-sqrt(var(y1));

emsigma2<-sqrt(var(y2));

mat<-matrix(6*51,51,6); # 결과물 저장 행렬

mat[1,1]<-0.5; #초기 값 저장

mat[1,2]<-0.5;

mat[1,3]<-emu1;

mat[1,4]<-emu2;

mat[1,5]<-emsigma1;

mat[1,6]<-emsigma2;

for (i in 1:50){ # do EM 50번 까지

w1<-ema1*dnorm(y, emu1, emsigma1);

w2<-ema2*dnorm(y, emu2, emsigma2);

nw1<-w1/(w1+w2); #normalize

nw2<-w2/(w1+w2);

ema1<-mean(nw1);

ema2<-mean(nw2);

emu1<-sum(nw1*y)/sum(nw1);

emu2<-sum(nw2*y)/sum(nw2);

emsigma1<-sqrt(sum(nw1*(y-emu1)^2)/sum(nw1));

emsigma2<-sqrt(sum(nw2*(y-emu2)^2)/sum(nw2));

mat[i+1,1]<-ema1; #EM 과정을 결과에 저장

mat[i+1,2]<-ema2;

mat[i+1,3]<-emu1;

mat[i+1,4]<-emu2;

mat[i+1,5]<-emsigma1;

mat[i+1,6]<-emsigma2;

}

davison모델선택책.pdf

4.73MB

davison모델선택책.pdf

4.73MB