Цикл решения ряда нелинейных уравнений с nleqslv в R

У меня проблема с решением ряда нелинейных уравнений с nleqslv в R, чтобы найти меру расстояния до значения по умолчанию. Это мой первый код R, поэтому я все еще борюсь с некоторыми проблемами. Мой код выглядит так (миниатюризированный до data.frame с тремя регистрами):

library("nleqslv")
D <- c(28000000, 59150000, 38357000)
VE <- c(4257875, 10522163.6, 31230643)
R  <- c(0.059883, 0.059883, 0.059883)
SE <- c(0.313887897, 0.449654737, 0.449734826976691)
df <- data.frame(D, VE, R, SE)

for(i in 1:nrow(df)){
  fnewton <- function(x){
    y <- numeric(2)
    d1 <- (log(x[1]/df$D[i])+(df$R[i]+x[2]^2/2))/x[2]
    d2 <- d1-x[2]
    y1 <- df$VE[i]-(x[1]*pnorm(d1)-exp(-df$R[i])*df$D[i]*pnorm(d2))
    y2 <- df$SE[i]*df$VE[i]-pnorm(d1)*x[2]*x[1]
    y
  }
  xstart <- c(df$VE[i], df$SE[i])
  df$VA[i] <- nleqslv(xstart, fnewton, method="Newton")$x[1]
  df$SA[i] <- nleqslv(xstart, fnewton, method="Newton")$x[2]
  i=i+1
}

Моя проблема в том, что мой код дает мне только одно решение, а это означает, что мой цикл не работает должным образом. Цикл должен преодолеть тот факт, что fnewton в первую очередь является вектором длины 2, но мои данные (или мой пример) являются более длинным вектором, чем 2. Я пробовал некоторые вещи, но не могу справиться с проблемой, я думаю, что есть простое решение для этого, но я не вижу своей ошибки.


r loops nonlinear-functions
person user78859    schedule 03.01.2016    source источник
comment
В функции fnewton должны ли y1,y2 быть y[1],y[2]?   -  person mra68    schedule 04.01.2016

Ответы (2)


arrow_upward
0
arrow_downward

В коде есть несколько простых ошибок.

1) Как прокомментировал mra68, измените y1, y2 в функции fnewton на y[1], y[2].

2) удалить последнюю строку i=i+1 (следующая i автоматически отображается строкой for(i in 1:nrow(df)){.)

3) (Необязательно) Инициализируйте df с указанными VA, SA.

Вот окончательный рабочий код:

library("nleqslv")
D <- c(28000000, 59150000, 38357000)
VE <- c(4257875, 10522163.6, 31230643)
R  <- c(0.059883, 0.059883, 0.059883)
SE <- c(0.313887897, 0.449654737, 0.449734826976691)
df <- data.frame(D, VE, R, SE, VA=numeric(3), SA=numeric(3))

for(i in 1:nrow(df)){

  fnewton <- function(x){
    d1 <- (log(x[1]/df$D[i])+(df$R[i]+x[2]^2/2))/x[2]
    d2 <- d1-x[2]

    y <- numeric(2)
    y[1] <- df$VE[i]-(x[1]*pnorm(d1)-exp(-df$R[i])*df$D[i]*pnorm(d2))
    y[2] <- df$SE[i]*df$VE[i]-pnorm(d1)*x[2]*x[1]
    y
  }

  xstart <- c(df$VE[i], df$SE[i])
  df$VA[i] <- nleqslv(xstart, fnewton, method="Newton")$x[1]
  df$SA[i] <- nleqslv(xstart, fnewton, method="Newton")$x[2]
}
person skwon    schedule 04.01.2016
comment
Неэффективность в этом коде при вызове nleqslv. Вам нужно только вызвать его один раз и присвоить результат переменной z, например. вот так z <- nleqlsv(....). А затем назначьте z$x[1] and z$x[2] соответствующим строкам df. Теперь вы дважды вызываете nleqslv с одинаковыми начальными значениями. - person Bhas; 05.01.2016

arrow_upward
0
arrow_downward

Код skwon можно сделать более эффективным, определив функцию fnewton вне цикла for и поместив переменные df и i в аргументы. Вот так

fnewton <- function(x,df,i){
  d1 <- (log(x[1]/df$D[i])+(df$R[i]+x[2]^2/2))/x[2]
  d2 <- d1-x[2]

  y <- numeric(2)
  y[1] <- df$VE[i]-(x[1]*pnorm(d1)-exp(-df$R[i])*df$D[i]*pnorm(d2))
  y[2] <- df$SE[i]*df$VE[i]-pnorm(d1)*x[2]*x[1]
  y
}

а затем измените цикл на этот

for(i in 1:nrow(df)){
    xstart <- c(df$VE[i], df$SE[i])
    z <- nleqslv(xstart, fnewton,  method="Newton",control=list(trace=1), df=df,i=i)
    df$VA[i] <- z$x[1]
    df$SA[i] <- z$x[2]
}

Если функция fnewton становится более сложной, вы также можете использовать пакет compiler, чтобы ускорить ее (немного).

Наконец, вы действительно должны проверить, что nleqslv действительно нашел решение, проверив, если z$termcd==1. Если нет, то пропустите присвоение значений z$x. Я оставлю это для вас.

person Bhas    schedule 05.01.2016