데이터 모형을 돌리고 나면, 모델에 대한 결과를 추출 해 올 때는 특별한 테크닉이 필요하다. 모형의 결과는 list화 되어 있기 때문에, 이를 어떻게 조정하고 추출 해오는지가 중요하다.
완벽하게 최적화 하기 위해서는 좌충우돌(左衝右突) 할 수 밖에 없다.
이번에는 은닉 마르코프 HMM 모형에서 데이터 최적화 하기 위한 도식이다.
 |
| HMM 은닉 마르코프 모델 데이터 서치 |
위 도식과 같이 소스코드를 정리한 것은 아래와 같다.
## 불필요 객체 지우기
rm(list = ls())
## 필요한 패키지 vector 형태 나열 및 함수형태로 불러오기
pkgs <- c("depmixS4", "ggplot2", "dplyr", "tqk", "forecast","stats")
sapply(pkgs, require, character.only = T)# 데이터가 바뀌지 않게 처리
set.seed(1004)
## 주식 가져오는 함수
data_get_fun <- function(days, stock_code){
# Data get from tqk
stock <- tqk_get(stock_code, from= days)
stock$종목코드 <- stock_code
# R은 데이터를 가져와 전저리 할때 단순함
stock <- stock %>%
arrange(date) %>%
mutate(전일비 = c(NA, diff(adjusted))) %>%
filter(!is.na(전일비)) %>%
mutate(전일시세 = adjusted - 전일비) %>%
mutate(등락율 = round((adjusted - 전일시세)/전일시세 * 100, 2))
return(stock)
}
## 데이터를 불러오고 HMM 모델링하는 함수
model_fit_fun <- function(state_count){
# HMM 모델 돌리기
hmm.model <- depmix(adjusted ~1, family = gaussian(),
nstates = state_count, data = df)
model.fit <- fit(hmm.model)
# 모델 결과 보기
fit <- attr(model.fit, "response")
# state 상태에 대한 조건부 확률 큰것
post_probs <- posterior(model.fit, type= "viterbi")
post_max <- as.matrix(post_probs[2:(state_count+1)])
post_prob_vt <- apply(post_max,1, max)
# 그래프를 그리기 위한 데이터 프레임 작업
df1 <- data.frame(date = df$date,
adjusted = df$adjusted,
state = post_probs$state,
post_prob = round(post_prob_vt,2))
# 확률이 95% 이하인것은 데이터 전처리
df1 <- df1 %>%
mutate(adjusted1 = ifelse(post_prob >=0.95, adjusted, NA ))
return(list(df1, fit))
}
## fitting 화 된 모델 데이터 프레임 만드는 function
state_fun <- function(i) {
state <- i
coeff <-as.double(fit[[i]][[1]]@parameters[["coefficients"]])
sd <- as.double(fit[[i]][[1]]@parameters[["sd"]])
df1 <- data.frame(state,
coeff = round(coeff,2),
sd = round(sd,2))
return(df1)
}
위에는 함수가 크게 3개가 있다. 주식 데이터를 가져오는 함수이고, 하나는 모델을 불러오는 함수이고, 모델을 데이터 프레임으로 변경 하는 함수이다. 코드 치는 방법이 여러가지가 있지만, 함수를 만드는 이유는 나중에 for 문이나 apply 계열의 패키지를 이용하여, 지속적으로 재사용 하려고 하기 때문이다.
우리가 원하는 데이터셋을 구하는 것은 state 와 확률이 들어간 Raw데이터 셋과 모델결과를 데이터 셋으로 만드는 것이다.
앞 편에서 주어진 코드는 데이터셋 모델 결과를 data.frame으로 구현하는 것처럼 보였는데, 문제가 생겨서 다시 아래와 같이 조정 하여 코드를 작성 하려 한다.
그리서 위의 함수를 먼저 만든 다음 나는 파레메터를 정한 다음 차례대로 가져오려고 한다.
## 기초변수 Parameter
# 분석 시작하고 싶은 일자
start_days = as.Date("2020-06-01")
# 증권종목코드(대상으로 준비함)
stock_code <- "001680"
# #데이터 가져오기
df <- data_get_fun(start_days,stock_code)
state_count <- 7
일단 파라메터를 정할 때, 분석 시작하고 싶은 일자와 종목코드, state 갯수 등을 정한 다음에, 아래와 같이 데이터를 돌리면, 데이터 프레임을 만들 수 있다.
## 그래프 분석할 데이터 셋 및 모델 가져오기
list <- model_fit_fun(state_count) df <- list[[1]]
fit <- list[[2]]
# state_df 가져오기 list 모형을 바꾸기 위해서는 Transpose 해야한다.
i <- 1:state_count
state_list <- t(sapply(i,state_fun))
state_df <- data.frame ( state = unlist(state_list[,1]),
coeff = unlist(state_list[,2]),
sd = unlist(state_list[,3]))
# state_prob 참일 가능성 찾기
state_prob <- df %>%
group_by(state) %>%
summarise( n = n(),
percent = table(is.na(adjusted1))["TRUE"] / n)
# state_df 다시 정리하기
state_df <- state_df %>%
left_join(state_prob, by ="state")
state_df ## state coeff sd n percent
## 1 1 26474.35 148.07 37 0.78378378
## 2 2 28046.02 180.32 26 0.53846154
## 3 3 25023.49 466.08 122 0.04098361
## 4 4 27628.46 140.68 31 0.45161290
## 5 5 27006.03 205.74 44 0.36363636
## 6 6 29000.48 586.34 31 0.19354839
## 7 7 25959.69 214.39 45 0.53333333모형을 만들게 되면, 기본적으로 리스트에서 데이터를 서브 세팅해서 가지고 오는데, 이것을 데이터 프레임으로 리턴 하더라도, sapply를 돌리게 되면, 다시 데이터 형태가 리스트 형태로 바뀐다.
이를 다시 풀 때는 리스트를 풀어서 다시 데이터 프레임을 만들어야 한다. 항상 sapply를 이용해서 연산 할 때는 항상 그러하다.
state별 기울기, 표준편차, state별 갯수, 그리고, state 별 조건부 확률이 95%이상인 경우의 통계량을 위와 같이 구할 수 있다.
다음 번에는 이것을 조금 더 간소화 할 수 있는 코드를 다시 찾아봐야 할 것 같다.
댓글 없음:
댓글 쓰기