Получить непрерывную модель как функцию двух независимых переменных вне их определенного диапазона - PullRequest
Новая
       95

Получить непрерывную модель как функцию двух независимых переменных вне их определенного диапазона

1 голос
/ 07 апреля 2020

У меня есть массив времен times=[58418 58422 58424 58426 58428 58430 58540 58542 58650 58654 58656 58658 58660 58662 58664 58666 58668 58670 58672 58674 58768 58770 58772 58774 58776 58778 58780 58782 58784 58786 58788 58790 58792 58794 58883 58884 58886 58888 58890 58890 58892 58894 58896 58898 58904]

и соответствующих наблюдаемых величин y_obs =[ 0.00393986 0.00522288 0.00820794 0.01102782 0.00411525 0.00297762 0.00463183 0.00602662 0.0114886 0.00176694 0.01241464 0.01316199 0.01108201 0.01056611 0.0107585 0.00723887 0.0082614 0.01239229 0.00148118 0.00407329 0.00626722 0.01026926 0.01408419 0.02638901 0.02284189 0.02142943 0.02274845 0.01315814 0.01155898 0.00985705 0.00476936 0.00130343 0.00350376 0.00463576 -0.00610933 0.00286234 0.00845177 0.00849791 0.0151215 0.0151215 0.00967625 0.00802465 0.00291534 0.00819779 0.00366089] и относительных ошибок:

y_obs_err = [6.12189334e-05 6.07487598e-05 4.66365096e-05 4.48781264e-05 5.55250430e-05 6.18699105e-05 6.35339947e-05 6.21108524e-05 5.55636135e-05 7.66087180e-05 4.34256323e-05 3.61131000e-05 3.30783270e-05 2.41312040e-05 2.85080015e-05 2.96644612e-05 4.58662869e-05 5.19419065e-05 6.00479888e-05 6.62586953e-05 3.64830945e-05 2.58120956e-05 1.83249104e-05 1.59433858e-05 1.33375408e-05 1.29714326e-05 1.26025166e-05 1.47293107e-05 2.17933175e-05 2.21611713e-05 2.42946630e-05 3.61296843e-05 4.23009806e-05 7.23405476e-05 5.59390368e-05 4.68144974e-05 3.44773949e-05 2.32907036e-05 2.23491451e-05 2.23491451e-05 2.92956472e-05 3.28665479e-05 4.41214301e-05 4.88142073e-05 7.19116984e-05]

Я определяю функцию, которая вычисляет значение y как функцию времени, количества параметров и другой независимой переменной

p= [ 2.82890497 3.75014266 5.89347542 7.91821558 2.95484056 2.13799544 3.32575733 4.32724456 8.2490644 1.26870083 8.91397925 9.45059128 7.95712563 7.58669608 7.72483557 5.19766853 5.93186433 8.89793105 1.06351782 2.92471065 4.49999613 7.37354766 10.11275281 18.94787684 16.40097363 15.38679306 16.33387783 9.44782842 8.29959664 7.07757293 3.42450524 0.93588962 2.515773 3.32857547 7.180216 2.05522399 6.06855409 6.1016838 10.8575614 10.8575614 6.94775991 5.76187014 2.09327787 5.88619335 2.62859611]

Здесь я определяю функцию для вычисления переменной y : 

import numpy as np
import matplotlib.pyplot as plt

from lmfit import Model
import scipy.integrate as it
import scipy.constants as scc
def new_f_function(t, sum, f0, a, b,c, T0):
    n= 2 * np.pi * sum 

    obs_f = f0 + it.cumtrapz(-a * p**c + b, t-T0, initial=0)
    new_f = obs_f*(1+sum/scc.c)
    return new_f

Я создаю модель и инициализирую свои параметры первым предположением: 

# Create Model
model = Model(new_f_function, independent_vars=['t'])

# Initialize Parameter
params = model.make_params()

params['sum'].value = 1.483 
params['sum'].min = 1.47
params['sum'].max = 1.50

.... # and so on for the others

, затем подгоняю модель 

`result = model.fit(y, params, weights=(1./y_err)), t=times, scale_covar=False)`

, чтобы получить наилучшие параметры result.best_fit. Наконец, я могу построить наилучшее соответствие с помощью

fig, (ax1, ax2) = plt.subplots(nrows=2, ncols=1, figsize=(3, 6), sharex='all') ax2.plot(times, result.best_fit, label='best fit')

Мой вопрос: как я могу построить / определить функцию в точках, где переменная p не определена (например, я нет данных)? Я предполагаю, что вопрос похож на этот: Scipy curve_fit: как построить подобранную кривую за пределами точек данных?

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

Ответы [ 3 ]

0 голосов
/ 07 апреля 2020

Ваш model объект хранит подходящую вам функцию, а наилучшие параметры для этого соответствия сохраняются в свойстве result.best_values.

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

new_times_to_evaluate = np.array([56789, 54321])  # Or whatever values you'd like

new_predictions = model.eval(  # evaluate your model
    **result.best_values,  # unpack the dictionary of best parameters to use for evaluation
    t=new_times_to_evaluate  # evaluate on these new times
)

ax2.plot(new_times_to_evaluate, new_predictions, label='my new values')
0 голосов
/ 08 апреля 2020

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

Чтобы сделать это, просто выполните: 

 result = model.fit(y, params, weights=(1./y_err)), t=times, 
                    scale_covar=False)
 predicted = result.eval(t=new_times)

Просто прочитайте это как «оцените результат в новое время».

0 голосов
/ 07 апреля 2020

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

Если у вас есть функция f (x, y) = z, у вас есть две независимые входные переменные x и y и одна выходная переменная (это ваша наблюдаемая переменная, если я прав).

Вы можете построить это, создав 2-мерный массив (сколько угодно, независимо от вашего времени до подгонки), передав каждую точку в функцию и построив все вычисленные z. Это не может быть построено в нормальном xy графике, очевидно. Один из способов - это 3D-график или 2-D цветная карта.

Не уверен, правильно ли я вас понял.

...