Fitting fractions of multivariate polynomial

Question

How to fit the multivariate fractional polynomial of the following form:

Given a function: y = f(x,z), a function of two variables x and z. More specifically it is of the form:

y = (x^2 + x^3)/(z^2 + z^3)

Numerator is a polynomial of a 3rd degree of a predictor x, and the denominator is also a polynomial of a 3rd degree of some predictor z.

I would like to fit polynomials for each of the predictor x and z, i.e. I need to find the coefficients A,B,C,D:

y = (A*x^2 + B*x^3)/(C*z^2 + D*z^3)

Basically, y is a ratio of two polynomials of degree 3. How to fit such a function?

Sample of a data frame below. I cannot post full data frame as it has more than 1000 rows.

y = c(-4.10594369806545, -4.23691458506868, -4.24690667936422, -3.53677470254628, -4.30406509320417, -4.19442802077908, -4.66857169733859, -2.82271310942235, -4.19720194766181, 3.52164353473802, -4.3917019001973, -5.41654474791269, 2.87471821731616, -3.85922481986118, -4.25370811223789, -3.57887855889961, -5.33913936106829, -4.11775265312012, -2.89958841300109, -4.18661983833127)

x = c(8.06526520889773, 9.39897529082673,9.07348918922699,7.5522372875608, 9.17294998275762,5.77455154554441, 9.2005930205213, 8.07309119969315, 7.42177579364465,8.18896686364888, 8.07868822922987, 8.50956416425175,9.71269017726113, 7.98378106897745, 7.69893619981345, 8.49576524400262, 8.02224091680654,8.25400859056484, 7.58171964012531, 8.35655484545343)

z = c(2.56494935746154, 4.99043258677874, 4.43081679884331,3.66356164612965,4.53259949315326,1.79175946922805,4.23410650459726, 5.38449506278909,3.13549421592915,4.34380542185368, 3.43398720448515,2.77258872223978,6.94985645500077,3.97029191355212, 3.40119738166216,4.39444915467244,2.19722457733622,3.91202300542815,4.06044301054642, 3.87120101090789)

dat = data.frame(cbind(y=y,x=x,z=z))

UPDATE:

Call to nls:

nls(y~(a*(x**2) + b*(x**3))/(c*(z**2) + d*(z**3)), dat, start=list(a=1,b=1,c=1,d=1))

With my extended dataset of more than 10^3 rows, I see the error I have just posted. With the sample that I provided here, I also get gradient problem. — aza07, Sep 26 '16 at 12:22
Thanks for your attempt. I still get the parameter initialization error for y and z. — aza07, Sep 26 '16 at 12:33
Did you go through the iterations of parameters? How did you get C=10? — aza07, Sep 26 '16 at 12:51
I can, however, find the best model by swiping through the parameters space and calculate the R^2 and other model goodness metrics. Do you think this can be the way to find the correct parameters? — aza07, Sep 26 '16 at 12:55
I think it whould be better to use `optim()` additionaly. `abc.f <- function(par) sum(((par[1]*x**2 + par[2]*x**3) / (par[3]*z**2 + z**3) - y)^2); optim(par=c(1, 1, -7), fn = abc.f, control=list(maxit = 1e+4, reltol=1e-12))` (example) — cuttlefish44, Sep 26 '16 at 13:15

Drey · Answer 1 · 2016-09-26T14:22:37.883

This is a nice problem that you have here. At least I learned some things from it. However, I have a feeling that this question resloves more around the solution of this specific task (university assignment?) than to be a general question.

But let me share an approach to the solution: What we have here

can be simplyfied to

Solving to y^\theta becomes numerically more managable. As we may see (and after trying very hard and failing to solve a non linear problem) it is actually a division of two linear models. So one approach is to estimate coefficients of the two problems separately. That is we fix coefficients a and b to find c and d and afterwards use c and d to find a and b.

Following code solves for coefficients

First some data

library(dplyr)

sampleData <- data.frame(x = runif(100, -100, 100), z = runif(100, -100, 100)) %>%
  mutate(y = ( (-2 * x^2) + (5 * x^3) ) /  (-4 * z^2 + 6 * z^3)) %>%
  mutate(zxfactor = z^2/x^2,
         yy = y * zxfactor)

now we solve for yy. With some random starting values...

init2 <- structure(runif(4, -10, 10), names=c("A", "B", "C", "D"))
coefab <- init2[c("A", "B")]
coefcd <- init2[c("C", "D")]

... we need to fit a linear model for a and b by

and a linear model for c and d by

# don't use for loop but determine a terminal condition... but i'm too lazy :-)
for(i in 1:100) {
  # make linear prediction using coeff. c and d
  sampleData <- sampleData %>%
    mutate(yab = yy * (coefcd[1] + coefcd[2] * z))
  # and fit a model for a and b
  coefab <- coef(lm(yab ~ x, sampleData))
  # then make a linear prediction using coeff. a and b
  sampleData <- sampleData %>%
    mutate(ycd = (coefab[1] + coefab[2] * x) / yy)
  # and fit a model for c and d
  coefcd <- coef(lm(ycd ~ z, sampleData))
} # repeat until satisfied

coefab
coefcd

Are we happy with the coefficients we have found ? Lets check:

optimFun <- function(params, out, x, z) {
  res <- (params[1] + params[2]*x)/(params[3] + params[4]*z)
  return( sqrt(sum( (out - res)^2 )) )
}

optimFun(c(coefab, coefcd), x = sampleData$x, z = sampleData$z, out = sampleData$yy)

> optimFun(c(coefab, coefcd), x = sampleData$x, z = sampleData$z, out = sampleData$yy)
[1] 1.951043e-12

Indeed we are since the difference between the modelled function estimates and data yy (the scaled one) is close to zero. Different iterations result in different parameter estimates since the problem is overdetermined. (Maybe someone can explain it in more detail)

Comments:

The estimates are closer to zero than estimates determined by nls
If you want to use optim you can use the optimFun. And for faster convergence you may even define derivate function
Its very nice problem to show where general optimization may fail and that it is always worth it to think about the problem at hand.
Try lm(yyz ~ x, data = sampleData %>% mutate(yyz = yy*(-4 + 6*z))) which will return exact values for a = -2 and b = 5. Given any two parameters you can find a matching pair that minimizes the function.

Fitting fractions of multivariate polynomial

1 Answers1