2017-01-20 30 views
0

Я нахожусь на уровне 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') 

ответ

0

Для

So basically WhatI'm asking is: how do I define a fixed c parameter value for each data set? c1=0.4; c2=1;c3=5, for example

Почему не просто делать

fit_params.add('c1', value=0.4, vary=False) 
fit_params.add('c2', value=1.0, vary=False) 
fit_params.add('c3', value=5.0, vary=False) 

после цикла for i in range(df.shape[1]/2):, определяющего другие параметры?

Для

What I'm doing wrong in the objective definition since when I run the code it states: TypeError: can't multiply sequence by non-int of type 'float'

Не видя полный отслеживающий (который, вероятно, сказать вам строку кода, где это исключение происходит), я бы предположил, что где-то есть список вместо Numpy массива умножается на поплавок.

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

+0

Я обновил информацию. Можете ли вы мне помочь, пожалуйста .... – Merrick

+0

Обновление полностью не связано с исходным вопросом, и вы даже не сказали, помогал ли более ранний ответ этому вопросу. –

+0

Обновление полностью не связано с исходным вопросом, и вы даже не сказали, помогал ли более ранний ответ. У Lmfit есть список рассылки: пожалуйста, используйте это. Без отслеживания и вывода и чтения кода мы можем только догадываться об этой проблеме. В дальнейших сообщениях дайте минимальный пример, полный текстовый вывод и полную трассировку. Я не понимаю, что вы подразумеваете под «36 кривой подходит». Существует одна целевая функция и один вызов для минимизации(): вы делаете одно подходящее. Вы зацикливаете на xs.shape [1] - это * должно быть * количество наборов данных, предположительно 6. Нет? –

 Смежные вопросы

  • Нет связанных вопросов^_^