CODEWITHSUNDEEP
rnonlinear-optimization
dslater0
dslater0

Reputation: 67

3 parameter nonlinear equation fitting in R

I'm trying to find a method in R which allows the parameters (a, b, c) of the general equation a*b^x + c which provides the best fit to 3 constrained random coordinates/points (p1, p2, p3 - with coordinates x1/y1, x2/y2 and x3/y3 respectively).

The constraints to these coordinates are:

I want to find a method which is able to generate values for a, b, and c which produces a line close to p1, p2 and p3. This is simply using desmos (see here for an example - https://www.desmos.com/calculator/4lmgazmrko) but I haven't been able to find a solution in R. I've tried the following:

x <- c(0, 0.7, 0.9)
y <- c(0.9, 0.8, 0)
df_test <- as.data.frame(cbind(x, y))

predict_y_nonlinearly <- function(beta, x){
  beta[1]*(beta[2]^x) + beta[3]
}

a_nonlinearmodel <- function(beta, x, y){
  y_hat <- predict_y_nonlinearly(beta, x)
  sum((y-y_hat)^2)
}

beta <- optim(rnorm(3), a_nonlinearmodel, method = "SANN",
              y = df_test$y, x = df_test$x)$par

predict_y_nonlinearly(beta, df_test$x)

But the optimisation function appears to get stuck in local minima, and rarely produces the correct solution (even when a different 'method' setting is used). I'm aware of the nls function, but this requires starting values to be chosen - I'd prefer a method which does not require manual input at this stage (as the desmos method is able to achieve).

Thanks

Upvotes: 0

Views: 201

Answers (1)

Ben Bolker
Ben Bolker

Reputation: 226162

Given the two zero constraints, we can reduce this analytically to a one-parameter problem:

x1 = 0 → y1 = a + c → c = y1-a
y3 = 0 → 0 = a*b^x3 + (y1-a)
→ a*(b^x3 - 1) = -y1
→ a = y1/(1-b^x3)

So we have a one-parameter function that predicts y, incorporating the x1 = y3 = 0 constraints:

predfun <- function(b = 1, x, y)  {
  a <- y[1]/(1-b^x[3])
  c <- y[1] - a
  a*b^x +c
}

A sum-of-squares target function:

target <- Vectorize(function(b) sum((y - predfun(b, x, y))^2))

Visualize:

curve(target, from = -10000, to = 100000, log = "y")

curve showing divergence at negative values, sharp minimum around 6e+04

Now use optimize() for 1D optimization (we still need to specify a starting interval, although not a specific starting point).

optimize(target, c(-10000, 1000000))

Results:

$minimum
[1] 58928.93

$objective
[1] 2.066598e-20

Upvotes: 2

Related Questions