Why do I get different power law models between R and Python?

73 Views Asked by At

I've noticed that when I make a power law model in R, it almost always under-predicts compared to one made in Python. Here is an R example:

#make data
x <- seq(1,10,1)
y <- c(1,2.1,3.5, 4.5, 4.8, 6, 9, 14, 20, 24.75)

model.pl <- lm(log(y) ~ log(x))
x_axis <- seq(1,10,1)
lines(x_axis, exp(predict(model.pl, data.frame(x=x_axis))), col = 'red', lty = 2, lwd=2)
summary(model.pl)

print(exp(model.pl$coefficients[1])) # 0.8093913 
print(model.pl$coefficients[2]) # 1.325091 

And in Python using Scipy:

lst = list(range(1,11,1))
lst2 = [1,2.1,3.5, 4.5, 4.8, 6, 9, 14, 20, 24.75]

df = pd.DataFrame(list(zip(lst, lst2)),
               columns =['x', 'y'])

def objective(x, a, b):
    return a * x ** b

x, y = df["x"], df["y"]

popt, _ = curve_fit(objective, x, y)

a, b = popt
print(a) #0.10104993 
print(b) #2.38573123 

So I end up with quite different parameters, which translates to different models when plotted below. Is this to be expected? I assumed they would fit essentially identical models:

Python model (left), R model (right)

As per phipsgabler, I modified the call, but this does not seem to change the model all that much:

model.pl <- lm(log(y) ~ log(x)-1)
x_axis <- seq(1,10,1)
lines(x_axis, exp(predict(model.pl, data.frame(x=x_axis))), col = 'red', lty = 2, lwd=2)
summary(model.pl)
print(exp(model.pl$coefficients[1]))
print(model.pl$coefficients[2])
0

There are 0 best solutions below