Fit to a power law

Question

I often encounter data which I hypothesize to be from a shifted power law,

$ y(x) = A x^k + B$.

I have in mind samples from an unknown deterministic function here, but you can think about a probability distribution if you prefer.

What is the best way to fit such data using Python? The shift means fitting a straight line in log-log space doesn't work.

Ideally I would prefer something a bit more convenient than using least squares directly, but if that is the best way, it is what it is. In that case, is a particular algorithm especially useful?

Answer 1

I would just spell out y(x) in pytorch and use autograd with squared loss together with one of the built in gradient descent algorithms on A, B and k. Unfortunately I don’t know any clever tricks to do it in a very simple way, if that was what you’re after. I had a go at it, here's the code

import torch
from torch.nn.functional import mse_loss
from collections import namedtuple
from itertools import count
from random import gauss, seed
seed(49)

# Our model
def f(x, A, B, k):
    return A*(x**k)+B

# Generate sample data
A = gauss(0, 10)
B = gauss(0, 10)
k = gauss(0, 10)
x = torch.arange(1, 10, dtype=torch.float)
y = f(x, A, B, k)
print(f"A={A}, B={B}, k={k}")
print(x)
print(y)

# Define the parameters
Params = namedtuple("Params", "A B k")

params = Params(A=torch.tensor([1.0], requires_grad=True),
                B=torch.tensor([1.0], requires_grad=True),
                k=torch.tensor([1.0], requires_grad=True))

# Fit the parameters
optimizer = torch.optim.Adam(params, lr=1)

tol = 1e-6 # MSE Loss error tolerance
for i in count():
    optimizer.zero_grad()
    y_hat = f(x, params.A, params.B, params.k)
    # Take the logarithm for improved numerical stability
    L = mse_loss(torch.log(y_hat), torch.log(y))
    if L < tol:
        break
    if i % 1000 == 0:
        print(f"Loss ({i}): {L.item()}")
    L.backward()
    optimizer.step()

print(params.A.item(), A)
print(params.B.item(), B)
print(params.k.item(), k)

Output:

A=9.42764034916727, B=4.212889473371394, k=12.835248103593624
tensor([1., 2., 3., 4., 5., 6., 7., 8., 9.])
tensor([1.3641e+01, 6.8901e+04, 1.2542e+07, 5.0349e+08, 8.8279e+09, 9.1657e+10,
        6.6290e+11, 3.6795e+12, 1.6686e+13])
Loss (0): 421.6722717285156
Loss (1000): 0.0005562781007029116
Loss (2000): 9.59674798650667e-06
9.461631774902344 9.42764034916727
4.17124080657959 4.212889473371394
12.833184242248535 12.835248103593624
```