Я нахожусь на уровне beginner в python, и я сталкиваюсь с проблемой. Я хотел бы использовать модуль lmfit для подгонки функции к 3 наборам данных (x, y) из CSV-файла с некоторыми общими параметрами (a, b, d) и одним фиксированным индивидуальным параметром (c). Мои данные в этом формате: x1 y1 x2 y2 x3 y3 .. .. .. .. .. .. .. .. .. .. .. .. .. ... . .. .. .. .. .. .. .. ..Установка с Lmfit нескольких datsets с общими параметрами и одним фиксированным параметром для каждого набора данных
.. .. .. .. .. ..
Вот пример показательной функции, и моя попытка подгонки к каждому набору данных и построению графика. У меня возникают проблемы с определением глобального остатка, а также для определения фиксированного значения для каждого набора данных параметра c. Это то, что я получил до сих пор:
import numpy as np
import pandas as pd
from lmfit import minimize, Parameters, report_fit
import matplotlib.pyplot as plt
#Load data
df=pd.read_csv('mydata.csv')
df[['x1', 'y1', 'x2', 'y2', 'x3','y3']][:6]
# set up the data
xs = np.array(df[df.columns[0::2]])
ys = np.array(df[df.columns[1::2]])
#define function
def biexp(xs,a,b,c,d):
return a*np.exp(-xs*b*c)+(1-a)*np.exp(-xs*d*c)
#Define a function that calculates biexp for data set i
def biexp_dataset(params, i, xs):
a = params['a_%i' % (i+1)].value
b = params['b_%i' % (i+1)].value
c = params['c_%i' % (i+1)].value
d = params['d_%i' % (i+1)].value
return biexp(xs,a,b,c,d)
#Define the real function to minimize, which calculates the total residual for all fits, modeled by biexp function
"""Where I'm stuck"""
def objective(params, xs,ys):
nys, nxs= df.shape[1]/2
resid = 0.0*ys[:]
# make residual per data set
for i in range(nys):
resid = ys[:, i] - biexp_dataset(params,i,xs)
# now flatten this to a 1D array, as minimize() needs
return resid.flatten()
#Define parameters
fit_params = Parameters()
for i in range(df.shape[1]/2):
fit_params.add('a_%i' % (i+1), value=1, min=0, max=9,vary=True)
fit_params.add('b_%i' % (i+1), value=0.3, min=0.0, vary=True)
fit_params.add('d_%i' % (i+1), value=0.5, min=0.0, vary=True)
"""How can I define c for each data sets?"""
"""Iwant for x1,y1 => c1=2; x2,y2=>c2=0.4; x3,y3=>c3=5, for example"""
# constrain all values of a,b and d to have the same value
for i in (2, 3):
fit_params['a_%i' % i].expr='a_1'
fit_params['b_%i' % i].expr='b_1'
fit_params['d_%i' % i].expr='d_1'
# run the global fit to all the data sets
result=minimize(objective, fit_params, args=(xs, ys))
report_fit(result.params)
# plot the data sets and fits
plt.figure()
for i in range(df.shape[1]/2):
y_fit = biexp_dataset(fit_params, i, xs)
plt.plot(xs[:,i], ys[:, i], 'o', xs[:,i], y_fit, '-')
plt.show()
Так в основном WhatI'm запрашиваемая: как мне определить значение параметра фиксированного с для каждого набора данных? с1 = 0,4; c2 = 1; c3 = 5, например
То, что я делаю неправильно в объективном определении, поскольку, когда я запускаю указанный код: TypeError: не может умножить последовательность на non-int типа 'float'
Любая помощь будет оценена ....
UpdateТак что я сделал какой-то новый код, используя реальные функции, которые я хотел, чтобы соответствовать (и что я уже установлена, используя глобальную опцию вписываться в OriginLab). Итак, мои наборы данных имеют одинаковые D1, D2, tau1 и tau2. Однако модель фиксируется как новая переменная, которая является экспериментальным временем (t_exps). Теперь у меня нет ошибок в коде, но то, что я получаю, - это не 6 встроенных кривых для шести наборов данных, а 36 кривых. Я заблудился ... Что случилось с моим кодом:
import numpy as np
import pandas as pd
from lmfit import minimize, Parameters, report_fit
import matplotlib.pyplot as plt
#set up the data
df=pd.read_excel('BTMS_hex.xlsx')
df = df.dropna() #drop rows that have invalid values (extra line in data file)
# set up the data
xs = np.array(df[df.columns[0::2]]).astype(float)
ys = np.log(np.array(df[df.columns[1::2]]).astype(float))
t_exps = [0.04857, 0.0989, 0.14923, 0.1993, 0.49957, 0.99957]
#Create function
def karger(xs,D1,D2,tau1,tau2,texp):
term1=D1+D2+(1/(xs/texp))*((1/tau1)+(1/(tau2)))
term2=np.sqrt(((D1-D2+(1/(xs/texp))*((1/tau1)+(1/(tau2))))**2)+(4/(((xs/texp)**2)*tau1*tau2)))
DA=0.5*(term1-term2)
DB=0.5*(term1+term2)
P1=tau1/(tau1+tau2)
P2=1-P1
CB=(P1*D1+P2*D2-DA)/(DB-DA)
return np.log((1-CB)*np.exp(-xs*DA)+CB*np.exp(-xs*DB))
def karger_dataset(params, xs,texp):
D1 = params['D1'].value
D2 = params['D2'].value
tau1 = params['tau1'].value
tau2 = params['tau2'].value
return karger(xs, D1,D2,tau1,tau2,texp)
def objective(params, xs,ys):
resid = np.array([])
for i in range(xs.shape[1]):
y_pred = karger_dataset(params, xs[:,i],t_exps[i])
resid=np.concatenate((resid,(ys[:, i] - y_pred)))
return resid
#Define parameters
fit_params = Parameters()
fit_params.add('D1', value=2E-1,min=0, vary=False)
fit_params.add('D2', value=1E-2, min=0, vary=True)
fit_params.add('tau1', value=8.2, min=0, vary=True)
fit_params.add('tau2', value=0.5, min=0, vary=True)
# run the global fit to all the data sets
result=minimize(objective, fit_params, args=(xs, ys))
report_fit(result.params)
print result.residual
# plot the data sets and fits
plt.figure()
for i in range(df.shape[1]/2):
y_fit = karger_dataset(fit_params, xs,t_exps[i])
plt.plot(xs[:,i], ys[:, i], 'o', xs[:,i], y_fit, '-')
"plt.yscale('log')"
plt.axis([0,1600, -10,0.1])
plt.savefig('books_read.png')
Я обновил информацию. Можете ли вы мне помочь, пожалуйста .... – Merrick
Обновление полностью не связано с исходным вопросом, и вы даже не сказали, помогал ли более ранний ответ этому вопросу. –
Обновление полностью не связано с исходным вопросом, и вы даже не сказали, помогал ли более ранний ответ. У Lmfit есть список рассылки: пожалуйста, используйте это. Без отслеживания и вывода и чтения кода мы можем только догадываться об этой проблеме. В дальнейших сообщениях дайте минимальный пример, полный текстовый вывод и полную трассировку. Я не понимаю, что вы подразумеваете под «36 кривой подходит». Существует одна целевая функция и один вызов для минимизации(): вы делаете одно подходящее. Вы зацикливаете на xs.shape [1] - это * должно быть * количество наборов данных, предположительно 6. Нет? –