pythoncurve-fittinggaussian

Multiple gaussian fit issue

I have this code that aims to fit the data in here https://drive.google.com/file/d/1uBrHxuftALiQTcTeBl-s6AHGiC8DGBWV/view?usp=sharing.

Where the function read_file is used to extract the information in a proper way. I do not what i am doing wrong that the gaussian fit is not right as you can see in the image (Gaussian fit) , and the standard deviation for all the fitted gaussian is the same or it varies very little (1e-09). Can anyone help me?

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

file_path='/content/drive/MyDrive/ProjetoXRF_Manuel/April2024/Amostra4.mca' # You need to change this to your directory

def read_file(file_path):
    # Abrir el archivo
    with open(file_path, 'r', encoding='latin1') as f:
        data = []
        live_time = None  # Inicializar variable live_time
        real_time = None  # Inicializar variable real_time
        for line in f:
            if 'LIVE_TIME' in line:
                live_time = line.strip()  # Almacenar la línea que contiene LIVE_TIME
            elif 'REAL_TIME' in line:
                real_time = line.strip()  # Almacenar la línea que contiene REAL_TIME
            elif all(char.isdigit() or char in ".-+e" for char in line.strip()):
                # Analizar y almacenar datos numéricos como flotante único
                row = [float(x) for x in line.split()]
                data.extend(row)

    # Convertir 'data' a un arreglo de numpy
    return np.array(data), live_time, real_time

data, live_time, real_time = read_numeric_data(file_path)
a = -0.0076702251954372525
b = 0.03952691936704189
Data=data

peaks, _ = find_peaks(Data[:, 0] / maxim, height=0.1)

# Get peak values and indices
peak_values = Data[peaks] / maxim
peak_indices = peaks

# Definición de la función gaussiana
def gaus(x, a, x0, sigma):
    return a * np.exp(-(x - x0) ** 2 / (2 * sigma ** 2))
limite=(10-a)/b
print(limite)
# Datos
x_values = a + b * np.linspace(0, 270, 270)
y_values = np.squeeze(Data[:270]) / np.max(Data[:270])

# Parámetros iniciales
sigma = 2 # initial value
standar=[]
# Ajuste de curvas para diferentes conjuntos de datos
for i in range(len(peak_values)):
    mean = a + b * peak_indices[i] 
    amplitude = peak_values[i][0]
    
    # Ajustar la curva para cada conjunto de datos
    popt, pcov = curve_fit(gaus, x_values, y_values, p0=[amplitude, mean, 1])
    
    # Visualización del ajuste
    m,n,p= popt
    plt.plot(x_values, gaus(x_values,amplitude,mean,p), label='fit')

    print(mean, amplitude,"Desviación standar",p)

    standar.append(p)

plt.scatter(x_values, y_values, label='Data')
plt.xlim(4, 10)
plt.xlabel('X')
plt.ylabel('Y')
plt.legend()
plt.grid(True)
plt.show()
Solution

  • The two main issues were:

    enter image description here

    The modified code:

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

file_path='Amostra4.mca'

def read_file(file_path):
    # Abrir el archivo
    with open(file_path, 'r', encoding='latin1') as f:
        data = []
        live_time = None  # Inicializar variable live_time
        real_time = None  # Inicializar variable real_time
        for line in f:
            if 'LIVE_TIME' in line:
                live_time = line.strip()  # Almacenar la línea que contiene LIVE_TIME
            elif 'REAL_TIME' in line:
                real_time = line.strip()  # Almacenar la línea que contiene REAL_TIME
            elif all(char.isdigit() or char in ".-+e" for char in line.strip()):
                # Analizar y almacenar datos numéricos como flotante único
                row = [float(x) for x in line.split()]
                data.extend(row)

    # Convertir 'data' a un arreglo de numpy
    return np.array(data), live_time, real_time

def multi_gaussian(x, *params):
    eps = 1e-10
    y = np.zeros_like(x)
    for i in range(0, len(params), 3):
        amplitude, mean, sigma = params[i:i+3]
        y += amplitude * np.exp(-(x - mean)**2 / (2 * sigma**2 + eps))
    return y

#Load and prepare data
data, live_time, real_time = read_file(file_path)

a = -0.0076702251954372525
b = 0.03952691936704189
y_start = 150
y_end = 270
data_y = data[y_start:y_end] / data[y_start:y_end].max()
data_x = a + b * np.linspace(y_start, y_end, num=len(data_y))

#Some initial estimates that could help with curve_fit() convergence
peak_indices, _ = find_peaks(data_y / data_y.max(), height=0.1)
peak_amplitudes = data_y[peak_indices]
peak_means = data_x[peak_indices]
peak_sigmas = [0.1] * len(peak_indices)

params_init = list(zip(peak_amplitudes, peak_means, peak_sigmas))
params_init = np.concatenate(params_init)

#Fit
params, params_cov = curve_fit(multi_gaussian, data_x, data_y, p0=params_init)

#Get a high-resolution interpolation of the fit
x_fine = np.linspace(data_x.min(), data_x.max(), num=500)
y_fit = multi_gaussian(x_fine, *params)

#
# Plot
#
fig, ax = plt.subplots(figsize=(8, 2))

#Data
ax.scatter(data_x, data_y, label='data', color='black', marker='o', s=20)

#Fitted curve
ax.plot(x_fine, y_fit, label='fit', linewidth=1.5, color='tab:purple')

#Annotate the estimated locations
[ax.axvline(mean, color='darkslategray', linewidth=1, linestyle=':') for mean in params[1::3]]

#Formatting
ax.set(xlabel='X', ylabel='Y', title='Data and fit')
ax.spines[['right', 'top']].set_visible(False)
ax.legend()