Экспоненциальная подгонка данных (python)

Question

Привет Я пытаюсь подогнать свои данные с помощью либо полиномиальной, либо экспоненциальной функции, с которой я потерпел неудачу в обоих. Код, я использую выглядит следующим образом:

with open('argon.dat','r') as f: 
    argon=f.readlines() 

eng1 = np.array([float(argon[argon.index(i)].split('\n')[0].split(' ')[0])*1000 for i in argon]) 
II01 = np.array([1-math.exp(-float(argon[argon.index(i)].split('\n')[0].split(' ')[1])*(1.784e-3*6.35)) for i in argon]) 

with open('copper.dat','r') as f: 
    copper=f.readlines() 

eng2 = [float(copper[copper.index(i)].split('\n')[0].split(' ')[0])*1000 for i in copper] 
II02 = [math.exp(-float(copper[copper.index(i)].split('\n')[0].split(' ')[1])*(8.128e-2*8.96)) for i in copper] 

fig, ax1 = plt.subplots(figsize=(12,10)) 
ax2 = ax1.twinx() 
ax1.set_yscale('log') 
ax2.set_yscale('log') 

arg = ax2.plot(eng1, II01, 'b--', label='Argon gas absorption at STP (6.35 cm)') 
cop = ax1.plot(eng2, II02, 'r', label='Copper wall transp. (0.81 mm)') 
plot = arg+cop 

labs = [l.get_label() for l in plot] 
ax1.legend(plot,labs,loc='lower right', fontsize=14) 
ax1.set_ylim(1e-6,1) 
ax2.set_ylim(1e-6,1) 
ax1.set_xlim(0,160) 
ax1.set_ylabel(r'$\displaystyle I/I_0$', fontsize=18) 
ax2.set_ylabel(r'$\displaystyle 1-I/I_0$', fontsize=18) 
ax1.set_xlabel('Photon Energy [keV]', fontsize=18) 
plt.show() 

Что дает мне , что я хочу сделать, это вместо данных чертежа, как это подогнать их к экспоненциальной кривой и умножать эти кривые, чтобы в итоге эффективность детектора (Я пытался умножить элемент за элементом, но у меня недостаточно точек данных, чтобы иметь плавную кривую). Я попытался использовать polyfit, а также попытался определить экспоненциальную функцию, чтобы увидеть ее работу, но в итоге я получил строку в обоих случаях

#def func(x, a, c, d): 
# return a*np.exp(-c*x)+d 
#  
#popt, pcov = curve_fit(func, eng1, II01) 
#plt.plot(eng1, func(eng1, *popt), label="Fitted Curve") 

и

В случае данных потребности берется из http://www.nist.gov/pml/data/xraycoef/index.cfm Также подобные работы могут быть найдены в Рис.3: http://scitation.aip.org/content/aapt/journal/ajp/83/8/10.1119/1.4923022

Покой eddited после @ ответ Оливера:

Я умножение с использованием существовавших данных :

i = 0 
eff1 = [] 
while i < len(arg): 
    eff1.append(arg[i]*cop[i]) 
    i += 1 

То, что я в конечном итоге это (красный цвет: медь, пунктирная синий: аргон, синий: умножение) Это то, что я предполагаю, чтобы получить, но с помощью функций для Thi кривых (будет комментарий, сделанный при ответе @ oliver относительно того, что неправильно или неправильно понял)

Answer 1

Причина, по которой curvefit дает вам постоянную (плоскую линию), заключается в том, что вы передаете ему набор данных, который некоррелирован с использованием модели, которую вы определили!

Позвольте мне воссоздать вашу установку первого:

argon = np.genfromtxt('argon.dat') 
copper = np.genfromtxt('copper.dat') 

f1 = 1 - np.exp(-argon[:,1] * 1.784e-3 * 6.35) 
f2 = np.exp(-copper[:,1] * 8.128e-2 * 8.96) 

Теперь обратите внимание, что f1 основан на 2-й колонке данных в файле argon.dat. Это не имеет отношение к первой колонке, хотя ничто не мешает вам черчение модифицированной версии 2-го столбца против первого, конечно, и это то, что вы делали, когда вы участок:

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

plt.semilogy(copper[:,0]*1000, f2, 'r-') # <- f2 was not based on the first column of that file, but on the 2nd. Nothing stops you from plotting those together though... 
plt.semilogy(argon[:,0]*1000, f1, 'b--') 
plt.ylim(1e-6,1) 
plt.xlim(0, 160) 

def model(x, a, b, offset): 
    return a*np.exp(-b*x) + offset 

Примечания: в модели у вас был параметр, называемый b, который не использовался. Это всегда плохая идея перейти к подгоночному алгоритму. Избавиться от этого.

Теперь вот трюк: вы сделали f1 на основе 2-го столбца, используя экспоненциальную модель. Поэтому вы должны передать curve_fit 2-й столбец в качестве независимой переменной (которая помечена как xdata в function's doc-string), а затем f1 в качестве зависимой переменной. Например:

popt1, pcov = curve_fit(model, argon[:,1], f1) 
popt2, pcov = curve_fit(model, cupper[:,1], f2) 

И это будет работать отлично.

Теперь, когда вы хотите построить гладкую версию продукта из двух графиков, вы должны начать с обычного интервала в независимой переменной. Для вас это энергия фотонов. 2-й столбец обоих файлов данных зависит от этого: есть функция (одна для аргона, другая для меди), которая связывает μ/ρ с энергией фотона.Итак, если у вас много данных для энергии, и вам удалось получить эти функции, у вас будет много данных для μ/ρ. Поскольку эти функции неизвестны, самое лучшее, что я могу сделать, это просто интерполировать. Однако данные логарифмичны, поэтому требуется логарифмическая интерполяция, а не линейная по умолчанию.

Итак, продолжайте получать много данных для энергии фотонов. В наборе данных energypoints экспоненциально растет, так что вы можете создать приличный новый набор точек с помощью np.logspace:

indep_var = argon[:,0]*1000 
energy = np.logspace(np.log10(indep_var.min()), 
        np.log10(indep_var.max()), 
        512) # both argon and cupper have the same min and max listed in the "energy" column. 

Он работает в нашу пользу, что энергия в обоих наборах данных имеют одинаковый минимум и максимум. В противном случае вам пришлось бы уменьшить диапазон этого лог-пространства.

Далее мы (логарифмически) интерполировать отношение energy -> μ/ρ:

interpolated_mu_rho_argon = np.power(10, np.interp(np.log10(energy), np.log10(indep_var), np.log10(argon[:,1]))) # perform logarithmic interpolation 
interpolated_mu_rho_copper = np.power(10, np.interp(np.log10(energy), np.log10(copper[:,0]*1000), np.log10(copper[:,1]))) 

Вот визуальное представление о том, что только что было сделано:

f, ax = plt.subplots(1,2, sharex=True, sharey=True) 
ax[0].semilogy(energy, interpolated_mu_rho_argon, 'gs-', lw=1) 
ax[0].semilogy(indep_var, argon[:,1], 'bo--', lw=1, ms=10) 
ax[1].semilogy(energy, interpolated_mu_rho_copper, 'gs-', lw=1) 
ax[1].semilogy(copper[:,0]*1000, copper[:,1], 'bo--', lw=1, ms=10) 
ax[0].set_title('argon') 
ax[1].set_title('copper') 
ax[0].set_xlabel('energy (keV)') 
ax[0].set_ylabel(r'$\mu/\rho$ (cm²/g)') 

Оригинальный набор данных, отмеченные с синими точками, была тонко интерполирована.

Теперь, последние шаги становятся легкими. Поскольку параметры вашей модели, которая отображает μ/ρ в какой-либо экспоненциальный вариант (функции, которые я переименовал как f1 и f2), уже были найдены, их можно использовать для создания гладкой версии данных, которые присутствовали, а также произведение этих двух функций:

plt.figure() 
plt.semilogy(energy, model(interpolated_mu_rho_argon, *popt1), 'b-', lw=1) 
plt.semilogy(argon[:,0]*1000, f1, 'bo ') 

plt.semilogy(copper[:,0]*1000, f2, 'ro ',) 
plt.semilogy(energy, model(interpolated_mu_rho_copper, *popt2), 'r-', lw=1) # same remark here! 

argon_copper_prod = model(interpolated_mu_rho_argon, *popt1)*model(interpolated_mu_rho_copper, *popt2) 
plt.semilogy(energy, argon_copper_prod, 'g-') 

plt.ylim(1e-6,1) 
plt.xlim(0, 160) 
plt.xlabel('energy (keV)') 
plt.ylabel(r'$\mu/\rho$ (cm²/g)') 

И там вы идете. Подводя итог:

1. генерировать достаточное количество точек данных независимой переменной, чтобы получить гладкие результаты
2. интерполировать отношения photon energy -> μ/ρ
3. сопоставить функцию интерполируемой μ/ρ