Подгонка параметров набора данных из модели - PullRequest
Новая
       87

Подгонка параметров набора данных из модели

0 голосов
/ 08 мая 2020

Я пытаюсь подогнать параметры транзитной кривой блеска.

Я наблюдал данные кривой транзитного блеска, и я использую .py в python, который через 4 параметра (период, a (большая полуось), наклон, pl anet радиус) возвращает модель кривая транзитного блеска. Я хотел бы минимизировать остаток между этими двумя кривыми блеска. Вот что я пытаюсь сделать: Во-первых, оценить максимальное правдоподобие, используя method = "L-BFGS-B", а затем применить mcm c с помощью emcee для оценки неопределенностей.

Код: 

p = lmfit.Parameters()
p.add_many(('per', 2.), ('inc', 90.), ('a', 5.), ('rp', 0.1))

per_b = [1., 3.]
a_b = [4., 6.]
inc_b = [88., 90.]
rp_b = [0.1, 0.3]

bounds = [(per_b[0], per_b[1]), (inc_b[0], inc_b[1]), (a_b[0], a_b[1]), (rp_b[0], rp_b[1])]

def residual(p):
    v = p.valuesdict()

    eclipse.criarEclipse(v['per'], v['a'], v['inc'], v['rp'])
    lc0 = numpy.array(eclipse.getCurvaLuz())   (observed flux data)
    ts0 = numpy.array(eclipse.getTempoHoras())   (observed time data)

    c = numpy.linspace(min(time_phased[bb]),max(time_phased[bb]),len(time_phased[bb]),endpoint=True)
    nn = interpolate.interp1d(ts0,lc0)

    return nn(c) - smoothed_LC[bb] (residual between the model and the data)

Внутренний остаток по умолчанию (p) Я убеждаюсь, что оба наблюдаемых данных (time_phased [bb] и Smoothed_LC [bb]) имеют одинаковый размер модели. кривая транзитного блеска. Я хочу, чтобы он давал мне наиболее подходящие значения для параметров (v ['per'], v ['a'], v ['in c'], v ['rp']).

Мне нужна ваша помощь, и я ценю ваше время и внимание. С наилучшими пожеланиями, Юрий.

1 Ответ

0 голосов
/ 09 мая 2020

Ваш пример неполный, в нем много частичных понятий, а некоторые недействительны Python. Это немного затрудняет понимание вашего намерения. Если приведенного ниже ответа недостаточно, обновите свой вопрос полным примером.

Кажется довольно очевидным, что вы хотите смоделировать свои данные smoothed_LC[bb] (не уверен, что такое bb) с моделью для некоторых эффект затмения. Исходя из этого предположения, я бы рекомендовал использовать подход lmfit.Model. Начните с написания функции, моделирующей данные, чтобы вы могли проверить и построить модель. Я не совсем уверен, что понимаю все, что вы делаете, но эта функция модели может выглядеть так: 

import numpy
from scipy import interpolate
from lmfit import Model
# import eclipse from somewhere....

def eclipse_lc(c, per, a, inc, p): 
    eclipse.criarEclipse(per, a, inc, rp)
    lc0 = numpy.array(eclipse.getCurvaLuz())   # observed flux data
    ts0 = numpy.array(eclipse.getTempoHoras()) # observed time data
    return interpolate.interp1d(ts0,lc0)(c)

С помощью этой функции модели вы можете построить модель: 

lc_model = Model(eclipse_lc)

, а затем создайте параметры для своей модели. Это автоматически назовет их после имен аргументов вашей модельной функции. Здесь вы также можете указать им начальные значения: 

params = lc_model.make_params(per=2, inc=90, a=5, rp=0.1)

Вы хотели установить верхнюю и нижнюю границы этих параметров. Это делается путем установки параметров min и max, а не упорядоченного массива границ: 

params['per'].min = 1.0
params['per'].max = 3.0

и так далее. Но также: установка таких жестких границ - обычно плохая идея. Установите границы, чтобы избежать нефизических значений параметров или когда становится очевидным, что вам нужно их разместить.

Теперь вы можете подогнать свои данные с помощью этой модели. Что ж, сначала вам нужно получить данные, которые вы хотите моделировать. Это кажется менее ясным из вашего примера, но возможно: 

   c_data = numpy.linspace(min(time_phased[bb]), max(time_phased[bb]),
                           len(time_phased[bb]), endpoint=True)
   lc_data = smoothed_LC[bb]

Ну, а зачем вам это c_data? Почему бы просто не использовать time_phased как независимую переменную? В любом случае, теперь вы можете подогнать свои данные к вашей модели с вашими параметрами: 

   result = lc_model(lc_data, params, c=c_data)

На этом этапе вы можете распечатать отчет о результатах и ​​/ или просмотреть или получить наиболее подходящие массивы: 

   print(result.fit_report())

   for p in result.params.items(): print(p)

   import matplotlib.pyplot as plt
   plt.plot(c_data, lc_data, label='data')
   plt.plot(c_data. result.best_fit, label='fit')
   plt.legend()
   plt.show()

Надеюсь, это поможет ...

...