Итеративная оптимизация альтернативного семейства glm - PullRequest
Новая
       8

Итеративная оптимизация альтернативного семейства glm

0 голосов
/ 17 апреля 2019

Я настраиваю альтернативную функцию ответа на часто используемую экспоненциальную функцию в пуассоновых гловах, которая называется softplus и определяется как $ \ frac {1} {c} \ log (1+ \ exp (c \ eta)) $, где $ \ eta $ соответствует линейному предиктору $ X \ beta $

. Я уже справился с оптимизацией, установив для параметра $ c $ произвольные фиксированные значения и выполняя только поиск $ \ hat {\ beta} $.

НО теперь для следующего шага я также должен оптимизировать этот параметр $ c $ (итеративно переключаясь между обновленным $ \ beta $ и текущим $ c $).

Я пытался написать функцию log-lik, функцию Score, а затем настроить оптимизацию Ньютона-Рафсона (используя цикл while), но я не знаю, как отделить обновление c от внешнего шага иобновление \ бета во внутреннем шаге ..

Есть предложения?

# Response function:
sp <- function(eta, c = 1 ) {  
  return(log(1 + exp(abs(c * eta)))/ c) 
} 

# Log Likelihood
l.lpois <- function(par, y, X){
  beta <- par[1:(length(par)-1)]
  c <- par[length(par)]
  l <- rep(NA, times = length(y))
  for (i in 1:length(l)){
    l[i] <- y[i] * log(sp(X[i,]%*%beta, c)) - sp(X[i,]%*%beta, c) 
  }
  l <- sum(l)
  return(l)
}

# Score function
score <- function(y, X, par){
  beta <- par[1:(length(par)-1)]
  c <- par[length(par)]

  s <- matrix(rep(NA, times = length(y)*length(par)), ncol = length(y))
  for (i in 1:length(y)){
    s[,i] <- c(X[i,], 1) * (y[i] * plogis(c * X[i,]%*%beta) / sp(X[i,]%*%beta, c) -     plogis(c * X[i,]%*%beta))
  }
  score <- rep(NA, times = nrow(s))
  for (j in 1:length(score)){
    score[j] <- sum(s[j,])
  }
  return(score)
}

# Optimization function
opt <- function(y, X, b.start, eps=0.0001, maxiter = 1e5){
  beta <- b.start[1:(length(b.start)-1)]
  c <- b.start[length(b.start)]

  b.old <- b.start
  i <- 0
  conv <- FALSE

  while(conv == FALSE){ 

    eta <- X%*%b.old[1:(length(b.old)-1)]
    s <- score(y, X, b.old)
    h <- numDeriv::hessian(l.lpois,b.old,y=y,X=X)

    invh <- solve(h)

    # update 
    b.new <- b.old + invh %*% s                                                         

    i <- i + 1

    # Test 
    if(any(is.nan(b.new))){                                                             
      b.new <- b.old                                                                
      warning("convergence failed")
      break 
    } 

    # convergence reached?
    if(sqrt(sum((b.new - b.old)^2))/sqrt(sum(b.old^2)) < eps | i >= maxiter){ 
      conv <- TRUE
    }
    b.old <- b.new
  }
  eta <- X%*%b.new[1:(length(b.new)-1)]

  # covariance
  invh  <- solve(numDeriv::hessian(l.lpois,b.new,y=y,X=X)) 


  fitted <- sp(eta, b.new[length(b.new)])

  result <- list("coefficients" = c(beta = b.new),
                 "fitted.values" = fitted,
                 "covariance" = invh)
}

# Running fails ..
n <- 100
x <- runif(n, 0, 1)
Xdes <- cbind(1, x) 
eta <- 1 + 2 * x
y <- rpois(n, sp(eta, c = 1))



opt(y,Xdes,c(0,1,1))

1 Ответ

0 голосов
/ 17 апреля 2019

У вас есть 2 ошибки:

строка 25: 

(y[i] * plogis(c * X[i,]%*%beta) / sp(X[i,]%*%beta, c) - plogis(c * X[i,]%*%beta))

это возвращает matrix, поэтому вы должны преобразовать в numeric: 

as.numeric(y[i] * plogis(c * X[i,]%*%beta) / sp(X[i,]%*%beta, c) - plogis(c * X[i,]%*%beta))

строка 23: ) отсутствует:

у вас есть: 

s <- matrix(rep(NA, times = length(y)*length(par), ncol = length(y))

пока должно быть: 

s <- matrix(rep(NA, times = length(y)*length(par)), ncol = length(y))
...