Power law fit doesn't work in python: it's either way off or returns only the starting parameters

Question

I'm very, very confused. I'm trying to fit a power law to my data. I tried my code to random generated data and it works just fine (see figure) but when I'm trying with my data, it's way off. I try to help curve_fit by giving starting values for the fitting parameters, but in this case, it's only returning the starting parameter for the fit. This does not make sense. Could anyone help me please ?

Here is my code:

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

# Define the power law function
def power_law(x, factor, exponent):
    '''
    x: x axis data
    factor: y axis intersection
    exponent: slope
    '''
    return factor * x ** exponent


# Define the power-law function

# Generate synthetic data following a power-law distribution
np.random.seed(0)  # for reproducibility
x_data = np.linspace(1, 10, 50)  # example x values
y_data = power_law(x = x_data, factor = 0.2, exponent = -10) * (1 + np.random.normal(scale=0.1, size=len(x_data)))  # example y values with added noise

# Fit the power-law model to the data
params, covariance = curve_fit(power_law, x_data, y_data)

# Extract fitted parameters
fac_fit, exp_fit = params

# Plot the data and the fitted power-law curve
plt.figure()
plt.scatter(x_data, y_data, label = 'Data')
plt.plot(x_data, power_law(x_data, fac_fit, exp_fit), color='red', label='Fitted Power Law')
plt.xscale('log')
plt.yscale('log')
plt.xlabel('X')
plt.ylabel('Y')
plt.title('Fitting a Power Law to Data')
plt.legend()
plt.grid(True)
plt.show()

# Print the fitted parameters
print("Fitted Parameters:")
print("factor =", fac_fit)
print("exponent =", exp_fit)

# My data
x_data = freq_full_ordered_withM[mask]
y_data =  PSD_full_ordered_withM[mask]

# Filter out non-positive values from x_data and corresponding y_data
positive_mask = x_data > 0
x_data = x_data[positive_mask]
y_data = y_data[positive_mask]

# Fit the power-law model to the data
params, covariance = curve_fit(power_law, x_data, y_data)

# Extract fitted parameters
fac_fit, exp_fit = params

# Plot the data and the fitted power-law curve
plt.figure()
plt.scatter(x_data, y_data, label = 'Data')
plt.plot(x_data, power_law(x_data, fac_fit, exp_fit), color='red', label='Fitted Power Law')
plt.xscale('log')
plt.yscale('log')
plt.xlabel('X')
plt.ylabel('Y')
plt.title('Fitting a Power Law to Data')
plt.legend()
plt.grid(True)
plt.show()

UPDATE: I gave it some starting parameter and now it's fitting something but it's still really off.

initial_guess = [10**2, -2]  # Initial guess parameters (a, b)
bounds = ([10**0, -5], [10**5, 0])  # Lower and upper bounds for parameters (a, b)

# Fit the power-law model to the data
params, covariance = curve_fit(power_law, x_data, y_data, p0=initial_guess, bounds=bounds)

Answer 1

For something to fit a power law on that decadal range of y values I think you would be better trying to fit a line for log(y) against log(x).

The correspondence is y = A.xⁿ corresponds to log(y) = log(A) + n * log(x)

The following does a better job than the original on your synthetic data, but I don't have access to your real data to try it.

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

# Define the power law function
def power_law(x, factor, exponent):
    return factor * x ** exponent

def linear_law(x, c, m):
    return m * x + c

# Generate synthetic data
np.random.seed(0)  # for reproducibility
x_data = np.linspace(1, 10, 50)        # example x values. NOTE: not well-distributed on a log scale
y_data = power_law(x = x_data, factor = 0.2, exponent = -10) * (1 + np.random.normal(scale=0.1, size=len(x_data)))  # example y values with added noise

params, covariance = curve_fit(linear_law, np.log( x_data ), np.log( y_data ) )

# Convert parameters back: y = Ax^n   <--->   log(y) = log(A) + n * log( x )
fac_fit, exp_fit = np.exp( params[0] ), params[1]

# Plot the data and the fitted power-law curve
plt.figure()
plt.scatter(x_data, y_data, label = 'Data')
plt.plot(x_data, power_law(x_data, fac_fit, exp_fit), color='red', label='Fitted Power Law')
plt.xscale('log');   plt.yscale('log')
plt.xlabel('X')  ;   plt.ylabel('Y')
plt.title('Fitting a Power Law to Data')
plt.legend()
plt.grid(True)
plt.show()

# Print the fitted parameters
print("Fitted Parameters:")
print("factor =", fac_fit)
print("exponent =", exp_fit)

Output:

Fitted Parameters:
factor = 0.2249726345273465
exponent = -10.07098332868395