Python中高斯函数拟合原子光谱发射线数据的方法及问题排查

Hey there! Let's break down why your current fitting approach isn't working, and how to fix it properly for your emission line data.

Why Your Current Code Isn't Working #

The scipy.stats.norm.fit function you're using is designed to estimate the parameters of a normal probability distribution from a set of sample points. It only looks at your wavelength values (he3888_1[:,0]) as if they're random samples from a normal distribution—it completely ignores your intensity data (he3888_1[:,1]).

That long low-intensity interval you suspected? It's absolutely throwing off the fit: all those low-intensity wavelength points are being treated as equal-weight samples, pulling the estimated mean (peak wavelength) away from the actual emission peak, and making the standard deviation much larger than it should be. On top of that, norm.pdf outputs a probability density (which sums to ~1), so it will never match the scale of your intensity values, making the fit line look totally disconnected from your data.

The Correct Approach: Curve Fitting to (Wavelength, Intensity) Data #

Instead of treating your wavelengths as a sample distribution, you need to fit a Gaussian function directly to your (x, y) data (wavelength vs. intensity). scipy.optimize.curve_fit is perfect for this—it minimizes the difference between your observed intensity values and the Gaussian function's predictions, using your intensity data as the target.

Here's a step-by-step implementation:

1. Define the Gaussian Function

First, write a standard Gaussian function that takes wavelength (x), amplitude (peak intensity), mean (peak wavelength), and standard deviation (width of the peak) as inputs:

def gaussian(x, amp, mean, std):
    return amp * np.exp(-(x - mean)**2 / (2 * std**2))

2. Full Fitting Code

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

# Load your data (assuming he3888_1 is your (wavelength, intensity) array)
x_data = he3888_1[:, 0]
y_data = he3888_1[:, 1]

# Define Gaussian function
def gaussian(x, amp, mean, std):
    return amp * np.exp(-(x - mean)**2 / (2 * std**2))

# Make an initial guess for parameters (critical for good convergence!)
# - amp: peak intensity from your data
# - mean: wavelength where intensity is maximum
# - std: rough estimate of peak width (e.g., 1/10 of the total wavelength range)
initial_guess = [
    np.max(y_data),
    x_data[np.argmax(y_data)],
    (np.max(x_data) - np.min(x_data)) / 10
]

# Perform the curve fit
popt, pcov = curve_fit(gaussian, x_data, y_data, p0=initial_guess)

# Extract fitted parameters
fit_amp, fit_mean, fit_std = popt

# Generate the fitted curve
x_fit = np.linspace(np.min(x_data), np.max(x_data), 100)
y_fit = gaussian(x_fit, fit_amp, fit_mean, fit_std)

# Plot results
plt.plot(x_data, y_data, color='r', label='Raw Emission Data')
plt.plot(x_fit, y_fit, color='b', linewidth=2, 
         label=f'Gaussian Fit\nPeak Wavelength: {fit_mean:.2f} Å')
plt.xlabel("Wavelength (Angstroms)")
plt.ylabel("Intensity")
plt.legend()
plt.show()

# Print key results
print(f"Fitted Peak Wavelength: {fit_mean:.2f} Angstroms")
print(f"Fitted Peak Intensity: {fit_amp:.2f}")
print(f"Peak Standard Deviation: {fit_std:.2f}")

Pro Tips for Better Fitting #

• Crop your data: If you have a lot of low-intensity background, crop the dataset to only include points around the emission peak (e.g., intensity > 10% of the maximum intensity). This removes noisy background that can skew the fit:

  # Crop to peak region
  peak_threshold = 0.1 * np.max(y_data)
  mask = y_data > peak_threshold
  x_cropped = x_data[mask]
  y_cropped = y_data[mask]
  
  # Use cropped data for fitting
  popt, pcov = curve_fit(gaussian, x_cropped, y_cropped, p0=initial_guess)
  

• Refine initial guesses: If your fit still doesn't converge, tweak the initial guess for std to be closer to the actual peak width (you can estimate this from your plot).

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