Как выполнить подгонку кривой семейства кривых с помощью python / scipy

Question

Мне нужно решить следующую проблему:

• Я хочу согласовать семейство кривых с одним уравнением (см. Код)
• Подогнанные параметры (C1-C4) должны быть постоянными
• Таким образом, я могу уместить (или описать) одну кривую семьи, просто изменив "Сигма и Т"

Я надеюсь, что смогу описать мою проблему, надеюсь, что вы, ребята, можете помогите, я был бы очень признателен!

Вопрос отредактирован (из-за недопонимания - 2020_04_04)

Сейчас я попытаюсь быть более точным c сейчас, к которому я прикрепил картинку, где Вы можете увидеть пример «семейства кривых», которое меняется для разных «сигм». Я хочу описать эти семейства кривых с помощью пары констант - C1, C2, C3 и C4 без их изменения. Подсказка состоит в том, чтобы найти оптимум констант, которые могут описать это семейство кривых, просто меняя сигма и Т в качестве переменных. Поэтому я должен подгонять параметры для множества кривых с минимумом ошибок. После этого уравнение должно охватить все семейство кривых, просто изменив «Сигма и Т».

Пример кривых

С наилучшими пожеланиями!

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


#Equation --> Eps_Cr = (C1*Sigma**C2*x**(C3+1)*e(-C4/T))/(C3+1)

def func(x, C1, C2, C3,C4):
    Sigma = 20
    T = 1
    return (C1*Sigma**C2*x**(C3+1)*np.exp(-C4*1/T))/(C3+1)

#Example Data 1
xdata = [1, 10, 100, 1000, 10000, 100000]
ydata = [0.000382,0.000407,0.000658,0.001169,0.002205,0.004304]

#Example Data 2
xdata1 = [1, 10, 100, 1000, 10000, 100000]
ydata1 = [0.002164,0.002371,0.004441,0.008571,0.016811,0.033261]

#Example Data 3
xdata2 = [1, 10, 100, 1000, 10000, 100000]
ydata2 = [0.001332,0.001457,0.002707,0.005157,0.010007,0.019597]

plt.plot(xdata, ydata, 'b-', label='data')
plt.plot(xdata1, ydata1, 'g-', label='data')
plt.plot(xdata2, ydata2, 'y-', label='data')

popt, pcov = curve_fit(func, xdata, ydata)

plt.plot(xdata, func(xdata, *popt), 'r--',
         label='fit: C1=%5.2e, C2=%5.3f, C3=%5.3f,C4=%5.3f' % tuple(popt))


plt.xlabel('X')
plt.ylabel('Y')
plt.legend()
plt.show()

Answer 1

Из дополнительной информации, которую вы предоставили в ответе, кажется, что вы хотите соответствовать иерархической модели. По крайней мере, так часто называют их статистики. Некоторые параметры совместно используются всеми точками данных (параметры от C1 до C4, а некоторые параметры совместно используются в группах наборов данных (T и Sigma). Все эти параметры необходимо оценивать по данным.

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

Поскольку вы уже используете curve_fit, я сделал версию вашего кода, которая выполняет свою работу. Стиль кода оставляет желать лучшего, так как я нет эксперта в scipy, но я думаю, что вы, по крайней мере, поймете метод.

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

def func(x_and_grp, C1, C2, C3, C4, Sigma0, Sigma1, Sigma2, T0, T1, T2):
    # We estimate one sigma and one T per group of data points
    x = x_and_grp[:,0]
    grp_id = x_and_grp[:,1]
    # here we select the appropriate T and Sigma for each data point based on their group id
    T = np.array([[T0, T1, T2][int(gid)] for gid in grp_id])
    Sigma = np.array([[Sigma0, Sigma1, Sigma2][int(gid)] for gid in grp_id])
    return (C1*Sigma**C2*x**(C3+1)*np.exp(-C4*1/T))/(C3+1)

#Example Data in 3 groups
xdata0 = [1, 10, 100, 1000, 10000, 100000]
ydata0 = [0.000382,0.000407,0.000658,0.001169,0.002205,0.004304]
xdata1 = [1, 10, 100, 1000, 10000, 100000]
ydata1 = [0.002164,0.002371,0.004441,0.008571,0.016811,0.033261]
xdata2 = [1, 10, 100, 1000, 10000, 100000]
ydata2 = [0.001332,0.001457,0.002707,0.005157,0.010007,0.019597]

#  merge all the data and add the group id to the x-data vectors
y_all = np.concatenate([ydata0, ydata1, ydata2])
x_and_grp_all = np.zeros(shape=(3 * 6, 2))
x_and_grp_all[:, 0] = np.concatenate([xdata0, xdata1, xdata2])
x_and_grp_all[0:6, 1] = 0
x_and_grp_all[6:12, 1] = 1
x_and_grp_all[12:18, 1] = 2

# fit a model to all the data together
popt, pcov = curve_fit(func, x_and_grp_all, y_all)

xspace = np.logspace(1,5)
plt.plot(xdata0, ydata0, 'b-', label='data')
plt.plot(xdata1, ydata1, 'g-', label='data')
plt.plot(xdata2, ydata2, 'y-', label='data')
for gid,color in zip([0,1,2],['r','k','purple']):
    T = popt[4+gid]
    Sigma = popt[7+gid]
    x_and_grp = np.column_stack([xspace,np.ones_like(xspace)*gid])
    plt.plot(xspace,
             func(x_and_grp, *popt),
             linestyle='dashed', color=color,
             label='fit: T=%5.2e, Sigma=%5.3f' % (T,Sigma))

plt.xlabel('X')
plt.ylabel('Y')
plt.title('fit: C1=%5.2e, C2=%5.3f, C3=%5.3f,C4=%5.3f' % tuple(popt[0:4]))
plt.legend()
plt.show()

Вывод выглядит так:

Наконец Я хочу добавить, что curve_fit не очень подходит для этой задачи, если у вас много разных групп. Рассмотрим другую библиотеку, которая может быть полезной. Статмодели могут быть возможны. Одна альтернатива - достичь scipy.optimize.minimze в вместо, так как это дает вам больше гибкости. Вам необходимо выполнить оценку доверительного интервала вручную, хотя ...

Я также хочу добавить, что описанный выше метод слишком сложен, если вы знаете T и Sigma для каждой группы данных. В этом случае мы добавляем соответствующее значение Sigma и T к x-вектору вместо идентификатора группы.

Answer 2

Согласно вашему запросу, я могу понять, что вам нужно подбирать одно уравнение для трех разных наборов данных отдельно. Итак, я обновил ваш код для того же самого, сохранив сигма и T одинаковыми. Пожалуйста, посмотрите и дайте мне знать.

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


#Equation --> Eps_Cr = (C1*Sigma**C2*x**(C3+1)*e(-C4/T))/(C3+1)

def func(x, C1, C2, C3,C4):
    Sigma = 20
    T = 1
    return (C1*Sigma**C2*x**(C3+1)*np.exp(-C4*1/T))/(C3+1)

#Example Data 1
xdata = [1, 10, 100, 1000, 10000, 100000]
ydata = [0.000382,0.000407,0.000658,0.001169,0.002205,0.004304]

#Example Data 2
xdata1 = [1, 10, 100, 1000, 10000, 100000]
ydata1 = [0.002164,0.002371,0.004441,0.008571,0.016811,0.033261]

#Example Data 3
xdata2 = [1, 10, 100, 1000, 10000, 100000]
ydata2 = [0.001332,0.001457,0.002707,0.005157,0.010007,0.019597]

plt.plot(xdata, ydata, 'b-', label='data 1')
plt.plot(xdata1, ydata1, 'g-', label='data 2')
plt.plot(xdata2, ydata2, 'y-', label='data 3')

popt, pcov = curve_fit(func, xdata, ydata)
popt1, pcov1 = curve_fit(func, xdata1, ydata1)
popt2, pcov2 = curve_fit(func, xdata2, ydata2)

plt.plot(xdata, func(xdata, *popt), 'r.',
         label='fit for Data 1: C1=%5.2e, C2=%5.3f, C3=%5.3f,C4=%5.3f' % tuple(popt))

plt.plot(xdata1, func(xdata1, *popt1), 'r+',
         label='fit for Data 2: C1=%5.2e, C2=%5.3f, C3=%5.3f,C4=%5.3f' % tuple(popt1))

plt.plot(xdata2, func(xdata2, *popt2), 'r--',
         label='fit for Data 3 : C1=%5.2e, C2=%5.3f, C3=%5.3f,C4=%5.3f' % tuple(popt2))


plt.xlabel('X')
plt.ylabel('Y')
plt.legend(loc='upper left',prop={'size': 8})
plt.show()