Loop for solving a number of non-linear equations

问题:

I have a problem in solving a number of non-linear equations with nleqslv in R in order to solve for a distance-to-default measure. This is my first R code, so I am still struggling with some problems. My code looks like this (miniaturized to a three-case-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
}

My problem is, that my code only gives me one solution, meaning that my loop does not work properly in the first place. The loop should overcome the fact, that fnewton is a vector of length 2 in the first place, but my data (or my example) is a longer vector than 2. I tried some things but I cannot handle the problem, I think there is a simple solution for this, but I do not see my mistake.

回答1:

There are some simple mistakes in code.

1) As mra68 commented, change y1, y2 in fnewton function to y[1], y[2].

2) remove last line i=i+1 (Next i is automatically mapped by the line for(i in 1:nrow(df)){.)

3) (Optional) Initialize df with VA, SA specified.

Here's final working code:

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]
}

回答2:

The code by skwon can be made more efficient by defining the fnewton function once outside the for loop and by putting the variables df and i in the arguments. Like so

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
}

and then change the loop to this

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]
}

If the function fnewton becomes more complicated you can also use package compiler to speed it up (a little bit).

Finally you really should test that nleqslv has in fact found a solution by testing if z$termcd==1. If not then skip the assigning of the z$x values. I'll leave that for you.