我有这个代码,旨在将数据放入这里DriveFilelink

read_file函数用于以 struct 化方式提取信息.然而,我在高斯匹配方面遇到了挑战,这从高斯fit image中观察到的差异中可以看出.这个问题似乎源于对某些参数施加的限制,例如固定所有高斯的平均值和四个幅度中的三个.尽管存在这些限制,但它们是必要的,因为它们是基于现有信息的.

有人知道如何解决这个问题吗?为了更贴合.

代码如下:

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import find_peaks
from scipy.optimize import curve_fit
from google.colab import drive
import os
# Mount Google Drive
drive.mount('/content/drive')


# Function to read numeric data from a file
def read_numeric_data(file_path):
    with open(file_path, 'r', encoding='latin1') as f:
        data = []
        live_time = None
        real_time = None
        for line in f:
            if 'LIVE_TIME' in line:
                live_time = line.strip()
            elif 'REAL_TIME' in line:
                real_time = line.strip()
            elif all(char.isdigit() or char in ".-+e" for char in line.strip()):
                row = [float(x) for x in line.split()]
                data.extend(row)
    return np.array(data), live_time, real_time



file_path = '/content/drive/MyDrive/ProjetoXRF_Manuel/April2024/20240402_In.mca'
data, _, _ = read_numeric_data(file_path)

a = -0.0188026396003431
b = 0.039549044037714
Data = data





# Función para convolucionar múltiples gaussianas
def multi_gaussian(x, *params):
    
    eps = 1e-10
    y = np.zeros_like(x)
    amplitude_relations = [1,0.5538, 0.1673 , 0.1673*0.5185 ]
    meanvalues= [24.210, 24.002 , 27.276 , 27.238]
    amplitude = params[0]

    for i in range(4):
        sigma = params[i * 3 + 2]  # Get the fitted standard deviation for this Gaussian
        
        y += amplitude*amplitude_relations[i]* np.exp(-(x - meanvalues[i])**2 / (2 * sigma**2 + eps))
    return y


sigma=[]
area=[]

# Función para graficar el espectro de energía convolucionado
def plot_convolved_spectrum(Data, a, b,i, ax=None):

    maxim = np.max(Data)
    workingarray = np.squeeze(Data)

    # Definir puntos máximos
    peaks, _ = find_peaks(workingarray / maxim, height=0.1)
    peak_values = workingarray[peaks] / maxim
    peak_indices = peaks

    # Calcular valores de energía correspondientes a los picos
    energy_spectrum = a + b * peak_indices

    # Preparar datos para convolución
    data = workingarray[:885] / maxim
    data_y = data / data.max()
    data_x = a + b * np.linspace(0, 885, num=len(data_y))

    # Ajustar múltiples gaussianas al espectro de energía
    peak_indices2, _ = find_peaks(data_y, height=0.1)
    peak_amplitudes = [1,0.5538, 0.1673 , 0.1673*0.5185 ]
    peak_means = [24.210, 24.002 , 27.276 , 27.238]
    peak_sigmas = [0.1] * 4

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

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

    # Obtener una interpolación de alta resolución del ajuste
    x_fine = np.linspace(data_x.min(), data_x.max(), num=20000)

    y_fit = multi_gaussian(x_fine, *params)

    #  Data Graphic
    ax.scatter(data_x, data_y, color='black', marker='o', s=20, label = 'Data' )
    y_data_error =np.sqrt(workingarray[:885])
    ax.plot(data_x, data_y + y_data_error/maxim, color='black',linestyle='--')
    ax.plot(data_x, data_y - y_data_error/maxim, color='black',linestyle='--')
    # Fit Graphic

    ax.plot(x_fine, y_fit, label="Fit", linewidth=1.5)



    # Extraer desviaciones estándar y amplitudes
    sigmas_array = params[2::3]
    # Calcular sigma promedio
    sigma.append(np.mean(sigmas_array))

    # Configuración del gráfico
    ax.set_xlabel('Energy (KeV)')
    ax.set_ylabel('Normalized Data')
    ax.legend()
    ax.set_title('Convolved Energy Spectrum')

    # Imprimir información
    print("Standard deviations:", sigmas_array)






fig, ax = plt.subplots()
plot_convolved_spectrum(Data, a, b,1,ax=ax)
ax.set_xlim(22, 28)
plt.show()

100 我try 过设置界限和细化初始猜测以改善高斯匹配.然而,需要注意的是,使用边界并不等同于修复参数.我的目标是仅维持两个参数的自由度:第一高斯的幅度和所有高斯的标准差.在探索这些调整时,我正在努力在约束和灵活性之间取得平衡,以实现所需的配合准确性.

Note This is the model: model

高斯函数的参数定义如下:A代表幅度,x_0代表平均值,西格玛代表标准差.在分析能谱的背景下,我拥有有关平均值和控制幅度关系的先验知识.因此,我的目标是根据该已知信息限制某些参数(例如A和平均值).具体来说,我打算仅适合第一个幅度,然后利用已知的关系,例如第二个峰值的A2 = 0.1673 * A1.并适应相应的标准偏差.

100

图中第一个峰值的明显奇异性可能表明单高斯.然而,事实并非如此.代表这个能量峰值的高斯实际上是由两个高斯加在一起组成的.这种复杂性之所以出现,是因为该水平的能量排放包括两个不同的值,24.002和24.210.由于实验的分辨率,我们无法仅通过看到情节来区分它们.

出于所有这些原因,我正在试图修复一些参数.

如果有人对这里的数据有问题,那就是:

extracted data.txt

推荐答案

这是一个MCVE,展示了如何从信号的重要部分回归可变数量的峰值.

import itertools
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import signal, stats, optimize

我们为您加载数据并进行后处理以获取物理系列:

data = pd.read_csv("20240402_In.mca", encoding="latin-1", names=["y"]).loc[12:2058]
data["y"] = pd.to_numeric(data["y"])
data["y"] /= data["y"].max()
data["x"] = -0.0188026396003431 + 0.039549044037714 * data.index

data = data.loc[(20. <= data["x"]) & (data["x"] <= 30.), :]

我们创建的模型接受可变数量的峰值(由find_peaks个解决方案驱动):

def peak(x, A, s, x0):
    law = stats.norm(loc=x0, scale=s)
    return A * law.pdf(x) / law.pdf(x0)

def model(x, *parameters):
    n = len(parameters) // 3
    y = np.zeros_like(x)
    for i in range(n):
        y += peak(x, *parameters[(i * 3):((i + 1) * 3)])
    return y

我们识别峰值:

indices, bases = signal.find_peaks(data.y, prominence=0.0125)
#(array([ 67, 118, 195, 210]),
# {'prominences': array([0.01987281, 0.99920509, 0.18918919, 0.03338633]),
#  'left_bases': array([  1,   1, 179, 205]),
#  'right_bases': array([ 76, 160, 219, 219])})

这是关键点,从识别开始,我们做出有根据的猜测:

x0s = data.x.values[indices]
As = bases["prominences"]
ss = (data.x.values[bases["right_bases"]] - data.x.values[bases["left_bases"]]) / 8.
p0 = list(itertools.chain(*zip(As, ss, x0s)))
#[0.01987281399046105,
# 0.37077228785356864,
# 22.682348638047493,
# 0.9992050874403816,
# 0.7860372502495658,
# 24.699349883970907,
# 0.1891891891891892,
# 0.19774522018856988,
# 27.744626274874886,
# 0.033386327503974564,
# 0.06921082706599968,
# 28.337861935440593]

现在我们根据您的数据集回归模型:

popt, pcov = optimize.curve_fit(model, data.x, data.y, p0=p0)
#array([1.03735804e-02, 9.61270732e-01, 2.29030214e+01, 9.92651381e-01,
#   1.59694755e-01, 2.46561911e+01, 1.85332645e-01, 1.21882422e-01,
#   2.77807838e+01, 3.67805911e-02, 1.21890416e-01, 2.83583849e+01])

我们现在可以估计模型:

yhat = model(data.x, *popt)

从图形上看,它会导致:

enter image description here

这对于无约束的配合来说是非常令人满意的,但肯定不能满足您想要实施的约束.

您可以在帖子中发布准备好编码为数组的数据吗?

Python相关问答推荐

Django序列化器没有验证或保存数据

带有pandas的分区列上的过滤器的多个条件read_parquet

Pandas read_jsonfuture 警告:解析字符串时,to_datetime与单位的行为已被反对

如何从格式为note:{neighbor:weight}的字典中构建networkx图?

如何使用矩阵在sklearn中同时对每个列执行matthews_corrcoef?

剧作家Python没有得到回应

优化在numpy数组中非零值周围创建缓冲区的函数的性能

pandas DataFrame GroupBy.diff函数的意外输出

如何比较numPy数组中的两个图像以获取它们不同的像素

如何在Python中将returns.context. DeliverresContext与Deliverc函数一起使用?

将两只Pandas rame乘以指数

迭代嵌套字典的值

如何并行化/加速并行numba代码?

如何指定列数据类型

如何使用Numpy. stracards重新编写滚动和?

幂集,其中每个元素可以是正或负""""

Flask运行时无法在Python中打印到控制台

BeautifulSoup-Screper有时运行得很好,很健壮--但有时它失败了::可能这里需要一些更多的异常处理?

GPT python SDK引入了大量开销/错误超时

如何使用matplotlib查看并列直方图