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