[R] logistic regression anlysis of mass data

이전 admission 데이터에 대해 glm 으로 구하는 것이 아닌 log.likelihode 와 optim으로 직접 구해보자.

setwd('/Users/bagjeong-yong/desktop/bigdataA/data')
admission <- read.csv('admission.csv',header=TRUE)
y <- as.numeric(admission[,1])
n <- length(y)
X<-as.matrix(cbind(rep(1,n),admission[,-1]))

log.L<-function(b){
  p.x<-exp(X%*%b) / (1 + exp(X%*%b))
  log.like <- sum(y*log(p.x) +(1-y)*log(1-p.x))
  return(-1*log.like)
}
log.L(c(0,0,0,0))
r <- optim(c(0,0,0,0),log.L)

 

해당 optim function 에서 hassian = TRUE 를 추가하면  log.likelihode에 이차도함수에 b^를 대입한 값을 준다.

이를 이용해 데이터를 나누고 합치자.

 

y<-as.numeric(admission[,1])
n<-length(y)
X<-as.matrix(cbind(rep(1,n),admission[,-1]))

log.like<-function(b){
p.X<-exp(X%*%b)/(1+exp(X%*%b))
-1*sum(y*log(p.X)+(1-y)*log(1-p.X))
}
log.like(c(0,0,0,0))
optim(rep(0,4),log.like,hessian=TRUE)
# hessian 이란 2차 편미분한 값을 모은 행렬
# argmin으로 편미분 하니까 log.like’’(b^) 의 결과

x1<- X[1:200,]
y1<- y[1:200]
x2<- X[-(1:200),]
y2<- y[-(1:200)]
# (-1:200) 이랑 -1:200 의미 다름. 후자는 -1부터 200까지.
log.like<-function(b){
p.x1<-exp(x1%*%b)/(1+exp(x1%*%b))
-1*sum(y1*log(p.x1)+(1-y1)*log(1-p.x1))
}

l1 <-optim(rep(0,4),log.like,hessian=TRUE)
log.like<-function(b){
p.x2<-exp(x2%*%b)/(1+exp(x2%*%b))
-1*sum(y2*log(p.x2)+(1-y2)*log(1-p.x2))
}
l2 <-optim(rep(0,4),log.like,hessian=TRUE)
A<-l1$hessian+l2$hessian
b<-l1$hessian%*%l1$par + l2$hessian%*%l2$par
solve(A) %*% b

  # $par argmin

또 만약 처음 데이터로 로지스틱 회귀를 진행하고  hassian 과 hassian*par 을 저장하고 있다가 이를 저장한 채로 

새로운 데이터가 들어올때 기존 정보와 새롭게 계산된 정보의 결과값을 누적시키면 online updating problem이 해결된다.

 

저번 주식데이터 예제에서 최고가를 통해 종가가 시가보다 큰지 작은지 예측하는 모델을 만들고 이를 1년씩 update 해보자. 

install.packages(“quantmod”)
library(quantmod)
start <- as.Date("2001-01-01")
end <- as.Date("2021-01-01")
getSymbols("AAPL",src="yahoo",from=start,to=end)
y <- as.numeric((AAPL$AAPL.Close-AAPL$AAPL.Open)>0)
n<-length(y)
x<-as.matrix(cbind(rep(1,n),as.numeric(AAPL$AAPL.High)))

log.like<-function(b){
p.x<-exp(x%*%b)/(1+exp(x%*%b))
-1*sum(y*log(p.x)+(1-y)*log(1-p.x))
}

l1 <-optim(rep(0,2),log.like,hessian=TRUE)
hat.beta.mle <- l1$par

A<-matrix(0,2,2)
b<-matrix(0,2,1)
start <- as.Date(c("2001-01-01",  "2006-01-01", "2011-01-01", "2016-01-01"))
end <- as.Date(c("2006-01-01", "2011-01-01", "2016-01-01", "2021-01-01"))
for(t in 1:4){
getSymbols("AAPL",src="yahoo",from=start[t],to=end[t])
y <- as.numeric((AAPL$AAPL.Close-AAPL$AAPL.Open)>0)
n<-length(y)
x<-as.matrix(cbind(rep(1,n),as.numeric(AAPL$AAPL.High)))

log.like<-function(b){
p.x<-exp(x%*%b)/(1+exp(x%*%b))
-1*sum(y*log(p.x)+(1-y)*log(1-p.x))
}
l <-optim(rep(0,2),log.like,hessian=TRUE)
A<-A+l$hessian
b<-b+l$hessian%*%l$par
hat.beta.online <- solve(A) %*% b