# How to model a power function ($y = ax^b$) with two exponents in two regions using R and nls

I have been unsuccessfully trying to model a relationship between two measured variables described by two different power functions, on either side of a threshold. My question is how to best estimate this relationship with a model. The aims of this model is to find the threshold $x_0$ and interpolate as exactly as possible the values of $y$ close to $x_0$. Estimating the relationship further from $x_0$ is less important. Based on theory I derived the following function: $$f(x) = \begin{cases} a(x_0-x)^b+c;& x\leq x_0\\ \frac{c-y_0}{d(x-x_0)+1}+y_0; & x > x_0. \end{cases}$$ where $a>0, b,c>0,d>0, \text{and } x_0>0$ are unknown parameters. This function behaves like a power function below $x_0$, and decreases roughly as $1/x$ to the asymptotic value of $y_0$ in the region above $x_0$. The function looks roughly like this:

So far I have tried using non-linear least squares to estimate this model, but I am getting the "singular gradient matrix at initial parameter estimates" error. My data and code in R is as follows:

x <- c( 0.33, 0.35, 0.39, 0.44, 0.48, 0.53, 0.57, 0.63, 0.74, 0.99, 1.12, 1.23, 1.37)
y <- c(72354.00, 23578.20, 1863.40, 743.80, 113.00, 9.80, 7.38, 5.30, 5.22, 5.03, 4.74, 4.53, 4.32)


and the code for the model I tried fitting is:

starting.values <- c(a = 8, b = 17, c = 8, d = 1,y_0 = 3, x_0 = .55)

model2 <- nls(y~ifelse(px < x_0,a*(x_0-px)^b+c,(c-y_0)/(d*(px-x_0)+1)+y_0),data = data.frame(x,y),
start = starting.values)


I have been trying a variety of likely starting parameters without success, is the problem too few data points, or is the model impossible to evaluate this way?

I have been successful, however, in modeling the relationship with a much simpler function: $$g(x) = \begin{cases} ax^b;& x\leq x_0\\ dx^c; & x > x_0. \end{cases}$$

where $d = ax_0^{b-c}$ to ensure continuity at the threshold. I did this by first fitting the data after applying the log-log transform, and than using the resulting values as inputs to the final model evaluated with nls. The code is as follows:

# transforming the data
x.log <- log(x)
y.log <- log(y)

# fitting a broken regression line to the log-log data
starting.values.log <- c(a = -5, b = 10, c = -.1, x_0 = .55)
model1 <- nls(y.log ~ifelse(x.log < log(x_0),b*x.log + a, c*x.log+b+log(x_0)*(a-c)),
data = data.frame(x.log,y.log),
start = starting.values.log)


below is the plot of the resulting model on the log-log plot:

and now I use the obtained parameters to fit the function:

# fitting the actual model using the parameters found previously
starting.values <- c(a = exp(-8.6238),b = -17.7984, c = -.4418, x_0 = 0.555)
model2 <- nls(y~ifelse(x < x_0,a*x^b,a*x_0^((b-c))*x^c),data = data.frame(x,y),
start = starting.values)


The parameters are:

Parameters:
Estimate Std. Error t value Pr(>|t|)
a    3.653e-05  1.130e-05   3.233   0.0103 *
b   -1.931e+01  2.804e-01 -68.857 1.45e-13 ***
c   -4.816e-01  8.832e+01  -0.005   0.9958
x_0  5.341e-01  1.335e+00   0.400   0.6984
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 388.5 on 9 degrees of freedom

Number of iterations to convergence: 18
Achieved convergence tolerance: 3.328e-07


The main problems with this method is that a) I am unable to determine the exponent of the power law when $x \rightarrow x_0$ from below and b) it does not allow for an asymptote other than 0 in $x \rightarrow \infty$. The first problem is much more important to me.

My question is how to best model the relationship between x and y, bering in mind that the theory supports the first function? Is the problem with trying to evaluate the first function too few data points or is it impossible to evaluate and I should try to use the second method instead? If so, than is there any way to obtain the exponent of the power law describing y when $x \rightarrow x_0$ from below?

• Your data really don't look much like the function you drew. If you take the logs of both variables it maybe looks vaguely like it, but only sort of. – Glen_b Feb 22 '15 at 4:48
• I think your problem is mainly with the curve in the left half. If you guess a few possible values for x0, you can just fit the left-half and kind of work out values for a b and c, but you only have 4 or 5 points to work with and the values are very unstable. – Glen_b Feb 22 '15 at 5:12
• @Glen_b Do I undertand correctly that my primary problem is too few data points in the $x<x_0$ region and not that I'm trying to estimate the model incorrectly? – Chris Novak Feb 22 '15 at 10:56
• I can't tell whether your code is right or not; I decided to start by trying to estimate the two components separately (after guessing at the changeover point). But too few points for the number of parameters isn't the only problem - the data are just not like your model. – Glen_b Feb 22 '15 at 11:33
• Worse, it seems like there's a ridge in parameter space, and so the fit is very sensitive to tiny changes. – Glen_b Feb 22 '15 at 11:39