Search code examples
rnonlinear-functions

Solving a system of 50 non-linear equations - error in fn(par,...) : argument is missing


I am trying to solve a system of 50 non-linear equations as follows:

Given a vector, y, which contains 49 different values, I'd like to transfer these 49 values into slightly different values in another vector, say, x such that:

log(x[1], y[1]) = n
...
log(x[49], y[49]) = n
x[1] + ... + x[49] = 1

For reasons of clarity, in the above equations y[i] is meant to be the base of the logarithm.

I've written the following code which, however, does not seem to work:

library(xlsx)
rm(list=ls())
setwd("C:/Users/.../folder")

my_data <- read.xlsx("samplefile.xlsx", 1)
y <- matrix(0:0, nrow=49,ncol=1)

for(i in 1:49) {
  if(my_data[i,1]!=0) {
    y[i,1] = 1/my_data[i,1]
    }
  }

for(i in 1:49) {
  fn <- function(x,n) {
  dummy1 <- log(x[i],y[i])-n
  dummy2 <- sum(x[1:49])-1

  return(c(dummy1,dummy2))
}
}

guess <- matrix(0.5:0.5, nrow = 50, ncol = 1)

nleqslv(guess,fn)

I'd expect that it solves for for x[i] and n. However, I get the following error message:

"Error in fn(par, ...) : argument "n" is missing, with no default"

Edit: Formatting


Solution

  • There are two ways of solving your problem.

    The first uses nleqslv for solving a system of 50 equations with 50 variables: x[1] to x[49] and n. The function to be solved by nleqslv must return a vector containing 50 elements; the system of equations must be square.

    The second simplifies your problem to a single equation in n.

    First solution:

    library(nleqslv)
    fn <- function(z,y) {
    
        x <- z[1:49]
        n <- z[50]
    
        f <- numeric(50)
    
        for( k in 1:49 ) f[k] <- log(x[k],y[k]) - n
        f[50] <- sum(x) - 1
        f
    }
    

    Generate some values for y and starting values for x and n and try to solve

    y <- rep(2,49)
    guess <- c( rep(.5,49), 1)
    
    res <- nleqslv(guess,fn,y=y)
    res
    

    Check result:

    # this should sum to 1
    sum(res$x[1:49])
    # value of n
    res$x[50]
    

    Second solution:

    Make use of the fact that log(a,b) = n implies a <- b^n. So log(x[k],y[k])=n is equivalent to x[k]=y[k]^n. So the x[1:49] can be computed immediately. We only need to determine n from the restriction that the elements of x sum to 1. This implies a simple function

    f2 <- function(n,y) {
        x <- y^n  # gives values for x
        sum(x) - 1
    }
    

    And now use uniroot like this assuming that the values for n will lie between -10 and 10 and that the signs of f2 at the endpoints have opposite sign.

    uniroot(f2,c(-10,10), y=y)
    

    You can run this code as given and check that nleqslv and uniroot give the same result for n.

    Whether this will work for your data is up to you find out.