Find exponent of the data-fit equation, R

Question

I would like to know how to adjust this equation/model below to my observed data (which are simple) so find the exponent p of this model

y~x^(-p)

My data is:

y=c(1.1178329,1.0871448,1.0897010,1.0759255,1.0535190,0.8725332)
x=c(6,5,4,3,2,1)

I tried the following model, but the values do not change and the iterations do not proceed.

library(minpack.lm)
mod <- nlsLM( y ~ x^(-p),
 start = c(p = 0.01) , 
trace = TRUE, lower=c(0.01) , upper=c(1))

iterations...

It.    0, RSS =  0.0671647, Par. =       0.01
It.    1, RSS =  0.0671647, Par. =       0.01

I thank everyone who can help me with this question.

Answer 1

There is no problem with this adjustment or code. See what happens when I change the initial kick to 0.9:

y=c(1.1178329,1.0871448,1.0897010,1.0759255,1.0535190,0.8725332)
x=c(6,5,4,3,2,1)

library(minpack.lm)
mod1 <- nlsLM( y ~ x^(-p),
  start = list(p = 0.9) , 
  trace = TRUE, lower=0.1, upper=1)
It.    0, RSS =    2.99354, Par. =        0.9
It.    1, RSS =   0.246237, Par. =        0.1
It.    2, RSS =   0.246237, Par. =        0.1    

More iterations occur as expected. If I keep the initial fixed value and reduce the lower bound of the search grid to 0.00001, see what happens:

mod2 <- nlsLM( y ~ x^(-p),
              start = list(p = 0.1) , 
              trace = TRUE, lower=0.00001, upper=1)
It.    0, RSS =   0.246237, Par. =        0.1
It.    1, RSS =  0.0544138, Par. =      1e-05
It.    2, RSS =  0.0544138, Par. =      1e-05

Apparently, the ideal value of p for this data set is quite close to 0. Note that including the RSS value (Residual Sum of Squares) reduces of mod1 for mod2, indicating that the second adjustment was better than the first.

In short, there’s nothing wrong with the code. It’s just converging too fast.

It may also be that the value of p is not positive. It may be that the value of p which best fits the data does not belong to the range determined by your initial fit, which is between 0.1 and 1. If this is true, the value of p will converge to something at the boundary of the interval defined in the function call. You came to think of this hypothesis?

See the code below, for example:

mod3 <- nlsLM( y ~ x^(-p),
              start = list(p = 0.1) , 
              trace = TRUE, lower=-1, upper=1)
It.    0, RSS =   0.246237, Par. =        0.1
It.    1, RSS =  0.0218977, Par. = -0.0816128
It.    2, RSS =  0.0166599, Par. = -0.0607616
It.    3, RSS =  0.0166587, Par. = -0.0604426
It.    4, RSS =  0.0166587, Par. = -0.0604428

mod3 resulted in the lowest RSS of all. Even, it has more expensive of adjustment, because it took longer to arrive in a response, with more interactions and a smoother convergence in the value of the parameter.

Answer 2

Either I’m completely wrong or lm enough to determine p.

y ~ x -p <=> y ~ e -plog(x) <=> log(y) ~ -plog(x)

Then we adjust this last linear model.

fit <- lm(log(y) ~ 0 + log(x))
coef(fit)

p <- -coef(fit)
p
     log(x) 
-0.06054118

plot(x, y, log = "xy")
abline(fit)

Answer 3

I would solve using the function optim which serves to do arbitrary optimizations, date a loss function with respect to some parameters.

Here is an example:

y=c(1.1178329,1.0871448,1.0897010,1.0759255,1.0535190,0.8725332)
x=c(6,5,4,3,2,1)

opt <- optim(runif(1), function(p){
  sum((y - x^(-p))^2)
}, method = "L-BFGS-B")

p <- opt$par

The function optim received three parameters:

• The initial value of p
• The loss function. In this case we are minimizing the sum of the squares of the residues of this model.
• The method of optimization.

The final value was -0.06044205. You can make a curve graph adjusted like this:

library(dplyr)
library(ggplot2)
data_frame(x = seq(1, 10, length.out = 100), y = x^(-p)) %>%
  ggplot(aes(x = x, y = y)) + 
  geom_line(colour = "blue")

The advantage of this approach is that you can use different loss functions and relationships between variables. Also easy to include regularisations etc.