Я пытаюсь установить профили линий, как обнаружено с помощью спектрографа на ПЗС. Для удобства рассмотрения я включил демонстрацию, которая, если ее решить, очень похожа на ту, которую я действительно хочу решить.
Я посмотрел на это: https://stats.stackexchange.com/info/46626/fitting-model-for-two-normal-distributions-in-pymc и различные другие вопросы и ответы, но они делают что-то принципиально иное, чем то, что я хочу делать.
import pymc as mc
import numpy as np
import pylab as pl
def GaussFunc(x, amplitude, centroid, sigma):
return amplitude * np.exp(-0.5 * ((x - centroid) / sigma)**2)
wavelength = np.arange(5000, 5050, 0.02)
# Profile 1
centroid_one = 5025.0
sigma_one = 2.2
height_one = 0.8
profile1 = GaussFunc(wavelength, height_one, centroid_one, sigma_one, )
# Profile 2
centroid_two = 5027.0
sigma_two = 1.2
height_two = 0.5
profile2 = GaussFunc(wavelength, height_two, centroid_two, sigma_two, )
# Measured values
noise = np.random.normal(0.0, 0.02, len(wavelength))
combined = profile1 + profile2 + noise
# If you want to plot what this looks like
pl.plot(wavelength, combined, label="Measured")
pl.plot(wavelength, profile1, color='red', linestyle='dashed', label="1")
pl.plot(wavelength, profile2, color='green', linestyle='dashed', label="2")
pl.title("Feature One and Two")
pl.legend()
Мой вопрос: Может ли PyMC (или какой-то вариант) дать мне среднее значение, амплитуду и сигма для двух компонентов, используемых выше?
Обратите внимание, что функции, которые я фактически буду вписывать в мою реальную проблему, НЕ являются гауссовыми. Поэтому, пожалуйста, предоставьте пример, используя общую функцию (например, GaussFunc в моем примере), а не "встроенный" pymc.Normal().
Кроме того, я понимаю, что выбор модели - это еще одна проблема: так что с текущим шумом 1 компонент (профиль) может быть все, что статистически оправдано. Но я хотел бы посмотреть, что будет лучшим решением для компонентов 1, 2, 3 и т.д.
Я также не согласен с идеей использования PyMC - если scikit-learn, astroML или какой-либо другой пакет кажутся идеальными, пожалуйста, дайте мне знать!
ИЗМЕНИТЬ:
Я потерпел неудачу несколькими способами, но одна из вещей, которые, я думаю, были на правильном пути, была следующей:
sigma_mc_one = mc.Uniform('sig', 0.01, 6.5)
height_mc_one = mc.Uniform('height', 0.1, 2.5)
centroid_mc_one = mc.Uniform('cen', 5015., 5040.)
Но я не мог создать mc.model, который сработал.