Подгонка данных с профилем Voigt в lmfit на python - огромные ошибки
Я пытаюсь сопоставить некоторые данные RIXS с профилями Voigt (lmfit в Python), и я определил профиль Voigt следующим образом:
def gfunction_norm(x, pos, gwid):
gauss= (1/(gwid*(np.sqrt(2*np.pi))))*(np.exp((-1.0/2)*((((x-pos)/gwid))**2)))
return (gauss-gauss.min())/(gauss.max()-gauss.min())
def lfunction_norm(x,pos,lwid):
lorentz=(0.15915*lwid)/((x-pos)**2+0.25*lwid**2)
return (lorentz-lorentz.min())/(lorentz.max()-lorentz.min())
def voigt(x, pos, gwid, lwid, int):
step=0.005
x2=np.arange(pos-7,pos+7+step,step)
voigt3=np.convolve(gfunction_norm(x2, pos, gwid), lfunction_norm(x2, pos, lwid), mode='same')
norm=(voigt3-voigt3.min())/(voigt3.max()-voigt3.min())
y=np.interp(energy, x2, norm)
return y * int
Я использовал это определение вместо популярного определения профиля Voigt в Python:
def voigt(x, alpha, cen, gamma):
sigma=alpha/np.sqrt(2*np.log(2))
return np.real(wofz((x-cen+1j*gamma)/sigma/np.sqrt(2)))/(sigma*2.51)
потому что это дает мне больше ясности в отношении интенсивности пиков и т. д.
Теперь у меня есть пара спектров с 9-10 пиками, и я пытаюсь сопоставить их все с профилями Фойгта (именно так, как я это определил).
Теперь у меня есть пара вопросов:
Как вы думаете, подходит ли мне мое определение Фойгта? Какие (не) преимущества я получаю, используя свертка вместо приближенного второго определения?
В результате моей подгонки иногда я получаю безумно большие стандартные отклонения. Например, это наиболее подходящие параметры для одного из пиков:
int8: 0.00986265 +/- 0.00113104 (11.47%) (init = 0.05)
pos8: -2.57960013 +/- 0.00790640 (0.31%) (init = -2.6)
gwid8: 0.06613237 +/- 0.02558441 (38.69%) (init = 0.1)
lwid8: 1.0909e-04 +/- 1.48706395 (1363160.91%) (init = 0.001)
(интенсивность, положение, ширина Гаусса и Лоренца соответственно). Означает ли этот вывод, что этот пик должен быть чисто гауссовым?
- Я заметил, что большие ошибки обычно случаются, когда параметр наилучшего соответствия очень мал. Это как-то связано с алгоритмом Левенберга-Марквардта, который по умолчанию используется в lmfit? Я должен отметить, что у меня иногда возникает та же проблема, даже когда я использую другое определение профиля Фойгта (и не только для лоренцевских ширин). Вот часть кода (это часть более крупного кода, и она находится в цикле for, что означает, что я использую одни и те же начальные значения для нескольких похожих спектров):
model = Model(final)
result = model.fit(spectra[:,nb_spectra], params, x=energy)
print(result.fit_report())
"final" - это сумма многих профилей Voigt, которые я определил ранее.
Спасибо!
Ответы 1
Это кажется дубликатом или продолжением Подгонка Lmfit вызывает огромные неопределенности - пожалуйста, используйте ТАК вопрос по теме.
Как вы думаете, подходит ли мне мое определение Фойгта? Какие (не) преимущества я получаю, используя свертка вместо приближенного второго определения?
Что заставляет вас говорить, что второе определение является приблизительным? В некотором смысле все вычисления с плавающей запятой являются приблизительными, но функция Фаддеевой из scipy.special.wofz является аналитическим решением для профиля Фойгта. Выполнение свертки самостоятельно, вероятно, будет немного медленнее, а также является приближением (на уровне машинной точности).
Поскольку вы используете Lmfit, я бы рекомендовал использовать его VoigtModel, который облегчит вам жизнь: он использует scipy.special.wofz и имена параметров, которые позволяют легко переключаться на другие профили (скажем, GaussianModel).
Вы не привели очень полный пример кода (для справки, минимальная рабочая версия фактического кода более или менее ожидается на SO и настоятельно рекомендуется), но это может выглядеть примерно так
from lmfit.models import VoigtModel
model = VoigtModel(prefix='p1_') + VoigtModel(prefix='p2_') + ...
В результате моей подгонки иногда я схожу с ума от большого стандарта отклонения. Например, это наиболее подходящие параметры для одного из пики:
int8: 0.00986265 +/- 0.00113104 (11.47%) (init = 0.05) pos8: -2.57960013 +/- 0.00790640 (0.31%) (init = -2.6) gwid8: 0.06613237 +/- 0.02558441 (38.69%) (init = 0.1) lwid8: 1.0909e-04 +/- 1.48706395 (1363160.91%) (init = 0.001)
(интенсивность, положение, ширина Гаусса и Лоренца соответственно). Означает ли этот вывод, что этот пик должен быть чисто гауссовым?
Во-первых, это не может быть «сумасшедшим» стандартным отклонением — оно как бы зависит от данных и остальной подгонки. Возможно, значение int8 очень, очень мало и сильно перекрывается с другими пиками — оно может сильно коррелировать с другими переменными. Но это вполне может означать, что пик больше похож на гауссовский.
Поскольку вы анализируете данные рассеяния рентгеновских лучей, использование функции Фойгта, вероятно, частично оправдано идеей (утверждением, предположением, ожиданием?), что отклик материала будет давать профиль Гаусса, в то время как приборы (включая источник рентгеновского излучения ) дало бы лоренцево уширение. Это говорит о том, что ширина Лоренца может быть одинаковой для различных пиков или, возможно, параметризована как простая функция длин волн падающего и рассеянного света или значений q. То есть вы могли бы (и, может быть, лучше) ограничить значения лоренцевской ширины (ваш lwid, я думаю, или gamma в lmfit.models.VoigtModel), чтобы все они были одинаковыми.
Спасибо за ваш ответ! Определенно кажется, что существует большая корреляция между гауссовой и лоренцевской шириной везде, где неопределенность лоренцевской ширины слишком велика. Кроме того, вы случайно не знаете, меняется ли разрешение в спектре потерь энергии для RIXS?