Блог

Я пытаюсь оценить некоторые параметры для функции на основе исторических данных (y). Это логистическая регрессия с асинхронным обновлением (дельта).

Это рассматриваемая система:

Где (для наших целей) X[t] = y[t], n2 формулируется как (1-n1), а phi1 и phi2 формулируются как phi1 deltaphi (как phi2> phi1), u[t] является неявным термином ошибки, зафиксированным y_hat.

Данные (значения журнала):

 y <- c(-0.083522212, -0.080744273, -0.098003453, -0.090994700, -0.105010991, -0.112070623, -0.115681762, -0.143194134, -0.146458642, -0.139691305, -0.128929970, -0.118047088, -0.095065509, -0.082399946, -0.100997872, -0.092699501, -0.082550517, -0.050470001, 0.030390893,  0.131429122,  0.180958369,  0.212374498,  0.223668346,  0.226461144,  0.209141361,  0.195377626,  0.178487458,
0.201981948,  0.233653604,  0.245221474,  0.227886405,  0.141274238,  0.046683795, -0.047819717, -0.112630561, -0.203788442, -0.238529171, -0.211924261, -0.233738086, -0.241872522, -0.238041656, -0.230753558, -0.242931741, -0.231894162, -0.243119603, -0.233052377, -0.230820606, -0.225594126, -0.232095554, -0.244800121, -0.252265025, -0.241778694, -0.227898251, -0.242656156, -0.229516117, -0.216082812, -0.220941314, -0.211800617, -0.183642284, -0.165779424, -0.159285263, -0.147407410, -0.138607996, -0.130455753, -0.094857132, -0.039392141, -0.003361144,  0.076381508,  0.101627405,  0.103042608,  0.096997308,  0.100308333, 0.098658702,  0.083952591,  0.077534743,  0.064491677,  0.056002466,  0.082643906,  0.080460147,  0.090688462)
  

Параметры для оценки:

 phi1, deltaphi, beta, delta
  

Другие материалы:

 T <- 80
y_hat <- mat.or.vec(nr = T, nc = 1)
R = 1.019656
alpha_bar = 0 
n1 = mat.or.vec(nr = T, nc = 1)
  

Функция:

 f <- function(t, phi1, deltaphi, beta, delta, R, alpha_bar, y) {
  n1[t] <<- (delta*n1[t-1]) (1-delta)*(1/(1 exp(-beta*((y[t-1] alpha_bar-R*y[t-2])*((-deltaphi))*y[t-3]))))
  y_hat = (((n1[t]*phi1 (1-n1[t])*(phi1 deltaphi))/((R alpha_bar)))*y[t-1])
return((y[t]-y_hat)^2) 
}
func <- function(par) sum(sapply(4:T,f, par[1],par[2], par[3],par[4], R, alpha_bar, y)) 

fit <- optimx(c(0.9, 1.05, 0.05, 0.6),
              method = "nlm",
              func,
              hessian = TRUE)
  

Я думаю, что проблема здесь в n1, я перепробовал кучу трюков, чтобы заставить его соответствовать, поскольку это рекурсивная переменная без начальной точки, следовательно, супероператор. Мне немного повезло, но функция всегда рушится на мне, поэтому я думаю, что функция может быть каким-то образом неправильно определена. Существуют некоторые ограничения на бета и дельта — бета должна быть положительно отличной от нуля, а дельта должна быть между 0 и 1.

Я могу получить некоторые ответы с помощью разных оптимизаторов (но не с nlm ); также действительно похоже, что ваша целевая функция может быть переопределена, что разные оптимизаторы будут обрабатывать по-разному… Я попытался немного упростить целевую функцию.

настройка констант

 maxT <- 80
R <- 1.019656
alpha_bar <- 0 
  

функция прогнозирования

Я использовал with(as.list(par) ...) , чтобы избежать необходимости распаковывать вектор параметров. Это немного усложняет отладку и заставляет нас использовать <<- внутри выражения …

 predfun <- function(par) {
    y_hat <- n1 <- numeric(maxT)  ## pre-allocate vectors
    with(as.list(par), {
        for (t in (4:maxT)) {
            n1[t] <<- (delta*n1[t-1]) (1-delta)*
                (1/(1 exp(-beta*((y[t-1] alpha_bar-R*y[t-2])*
                                 ((-deltaphi))*y[t-3]))))
            y_hat[t] <<- (((n1[t]*phi1 (1-n1[t])*(phi1 deltaphi))/
                           ((R alpha_bar)))*y[t-1])
        }
    })
    return(y_hat)
}
  

целевая функция

 func <- function(par) {
    return(sum((y[4:maxT]-predfun(par)[4:maxT])^2))
}
  

оптимизация

 p0 <- c(phi1=0.9, deltaphi=1.05, beta=0.05, delta=0.6)
func(p0)  ## check to make sure the function works for starting values
library(optimx)
fit <- optimx(par=p0,
              method = c("Nelder-Mead","BFGS","nlm"),
              func,
              hessian = TRUE)
  

Результаты

                   phi1 deltaphi     beta       delta         value fevals
Nelder-Mead -0.4967063 3.026097 6.227032 -0.16820934  7.546323e-02    265
BFGS        -0.4462906 2.903375 2.554713 -0.01802873  7.922255e-02    110
nlm                 NA       NA       NA          NA 8.988466e 307     NA
            gevals niter convcode  kkt1  kkt2 xtime
Nelder-Mead     NA    NA        0  TRUE  TRUE 0.143
BFGS           100    NA        1 FALSE FALSE 0.494
nlm             NA    NA     9999    NA    NA 0.001
  

nlm выдает все значения NA; BFGS не соответствует обоим критериям KKT; Nelder-Mead отображается нормально. Нелдер-Мид и BFGS дают умеренно разные оценки параметров, хотя оба дают, по-видимому, небольшой SSQ.

 cc <- coef(fit)
plot(y)
lines(predfun(cc["Nelder-Mead",]), col=2)
lines(predfun(cc["BFGS",]), col=4)
  

Подогнанные линии (почти) неразличимы.

преобразованные параметры

Если вы хотите подогнать параметры в преобразованном масштабе, вы должны предоставить начальные параметры вашей функции оптимизации в преобразованном масштабе, например

 p0 <- c(phi1=0.9, deltaphi=1.05, log_beta=log(0.05), logit_delta=qlogis(0.6))
  

Затем в вашей целевой функции преобразуйте их обратно перед их использованием:

 predfun <- function(par) {
    y_hat <- n1 <- numeric(maxT)  ## pre-allocate vectors
    with(as.list(par), {
       delta <- plogis(logit_delta)
       beta <- exp(log_beta)
       ### ... then the rest of your objective function ...