基于Scipy的半对数图曲线拟合与插值技术咨询

Hey there, let's fix this semilog fitting issue properly. The key problem with your current approach is that nonlinear curve fitting (like curve_fit for the exponential function) is sensitive to initial guesses, and cubic interpolation is meant to pass through all data points—not to capture the underlying trend. Here's a more reliable method tailored for semilog plots, plus how to handle extrapolation:

Why Your Current Methods Aren't Working #

• Cubic interpolation: This is an interpolation technique, not a fitting method. It will wiggle through every data point, which is great for filling gaps but terrible for identifying trends when data has noise (which yours does).
• curve_fit for exponential: Your initial guess x0 = [1e-6, 1e-6] is way off—looking at your data, as Mass500 increases, y500 generally decreases, so the slope a in 10^(a*x + b) should be negative. A positive initial guess leads the optimizer to converge on a bad local minimum.

The Better Approach: Linear Regression on Log-Transformed Data #

A straight line on a semilog y-axis means log10(y) is linear with x. Instead of fitting the exponential directly, we can:

1. Take the base-10 logarithm of your y500 values.
2. Perform linear regression on log10(y500) vs Mass500—this gives us the optimal least-squares fit, no finicky initial guesses needed.
3. Convert the linear fit back to the original exponential form for plotting.

Modified Code #

import numpy as np
import matplotlib.pyplot as plt

# Your original data
Mass500 = np.array([13.938, 13.816, 13.661, 13.683, 13.621, 13.547, 13.477, 13.492, 13.237, 13.232, 13.07, 13.048, 12.945, 12.861, 12.827, 12.577, 12.518])
y500 = np.array([7.65103978e-06, 4.79865790e-06, 2.08218909e-05, 4.98385924e-06, 5.63462673e-06, 2.90785458e-06, 2.21166794e-05, 1.34501705e-06, 6.26021870e-07, 6.62368879e-07, 6.46735547e-07, 3.68589447e-07, 3.86209019e-07, 5.61293275e-07, 2.41428755e-07, 9.62491134e-08, 2.36892162e-07])

# Step 1: Log-transform the y-values
y_log = np.log10(y500)

# Step 2: Linear regression (polyfit returns [slope, intercept])
a, b = np.polyfit(Mass500, y_log, deg=1)

# Step 3: Define the fit function for semilog plot
def semilog_line(x):
    return 10 ** (a * x + b)

# Generate points for plotting (including extrapolation range)
x_fit = np.linspace(min(Mass500) - 0.5, max(Mass500) + 0.5, 100)
y_fit = semilog_line(x_fit)

# Plotting
plt.figure(figsize=(10,6))
plt.semilogy(Mass500, y500, 'o', label='Raw Data')
plt.semilogy(x_fit, y_fit, 'r-', linewidth=2, label=f'Optimal Fit: $y = 10^{{{a:.2f}x + {b:.2f}}}$')

# Optional: Plot your old methods for comparison
from scipy.interpolate import interp1d
f2 = interp1d(Mass500, y500, kind='cubic')
plt.semilogy(Mass500, f2(Mass500), '--', label='Cubic Interpolation')

plt.xlabel('Mass500')
plt.ylabel('y500')
plt.legend(loc='best')
plt.grid(True, which="both", ls="-")
plt.show()

# Print the fit parameters
print(f"Fit parameters: slope (a) = {a:.4f}, intercept (b) = {b:.4f}")

Handling Extrapolation #

Extrapolation is straightforward with this method: just plug in any x value (even outside your original data range) into the semilog_line function. For example:

# Extrapolate to Mass500 = 14.0
extrapolated_y = semilog_line(14.0)
print(f"Extrapolated y at Mass500=14.0: {extrapolated_y:.2e}")

# Extrapolate to Mass500 = 12.0
extrapolated_y_12 = semilog_line(12.0)
print(f"Extrapolated y at Mass500=12.0: {extrapolated_y_12:.2e}")

Why This Works #

• Linear regression via np.polyfit uses ordinary least squares (OLS), which is the optimal method for fitting linear relationships—no guesswork required.
• By transforming to log-space, we turn the exponential relationship into a linear one, making the fit stable and robust even with small datasets.

内容的提问来源于stack exchange，提问作者bhjghjh