lmfit для построения параметров - PullRequest
Новая
       7

lmfit для построения параметров

0 голосов
/ 21 июня 2019

Я пытаюсь провести анализ чувствительности простой системы химических реакций.A -> B (со скоростью реакции k1) и A1 -> B (k2), BC (k3), C-> B (k4).Дело в том, чтобы построить ставки.Поэтому, если это небольшое значение k1, должно быть большое значение k2, и если k3 и k4 должны быть равны.Для k1 и k2 я получаю какую-то тенденцию, но k4 всегда одинаков.Что мне не хватает?Я изменяю параметры неправильно?

Моя попытка 

import matplotlib.pyplot as plt
import numpy as np
from lmfit import Parameters, report_fit, Minimizer, minimize, printfuncs, conf_interval, conf_interval2d
from scipy.integrate import odeint
from sys import exit
time = 10
Nt = 200
tt = np.linspace(0,time, Nt)   
p = Parameters()
# add with tuples: (NAME VALUE VARY MIN  MAX  EXPR  BRUTE_STEP)
p.add_many(('k1', 0.5, True, 0,1), ('k2', 0.5, True, 0,1), ('k3', 0.5, True, 0,1), ('k4', 0.5, True, 0,1))

time = 10
Nt = 11
tt = np.linspace(0,time, Nt)


def f(xs, t, ps):
    """Test"""
    try:
        k1 = ps['k1'].value
        k2 = ps['k2'].value
        k3 = ps['k3'].value
        k4 = ps['k4'].value
    except:
        k1, k2, k3, k4 = ps

    a, a1, b, c = xs
    da_dt = - k1*a
    da1_dt = - k2*a1
    db_dt = k1*a + k2*a1 + k4*c - k3*b
    dc_dt = k3*b - k4*c
    return da_dt, da1_dt, db_dt, dc_dt

def g(t, x0, ps):
    """
    Solution to the ODE x'(t) = f(t,x,k) with initial condition x(0) = x0
    """
    x = odeint(f, x0, t, args=(ps,))
    return x

def residual(ps, ts, data):
    x0 = np.array([1,0.5,0,0])
    model = g(ts, x0, ps)
    return (model)


#exit()
x0 = np.array([1,0.5,0,0])

k1, k2, k3, k4 = 1,1,1,1
true_params = np.array((k1,k2,k3,k4))
data = g(tt, x0, true_params)
data += ((np.random.normal(size=data.shape))) 

# create Minimizer
mini = Minimizer(residual, p,fcn_args=(tt,data))

# first solve with Nelder-Mead
out1 = mini.minimize(method='emcee')

out2 = mini.minimize(params=out1.params, method='Powell')

print(fit_report(out2))
print(report_fit(out1.params, min_correl=2))
##
ci, trace = conf_interval(mini, out2, sigmas = [0.01,0.02,0.3,0.4],
                                trace=True, verbose=True, maxiter = 200)


figure(1)
x, y, prob = trace['k1']['k1'], trace['k1']['k2'], trace['k1']['prob']
#x, y, prob = trace['k2']['k2'], trace['k2']['k1'], trace['k2']['prob']
scatter(x, y, c=prob ,s=30)
#plt.gca().set_xlim((0.00, 0.1))
#plt.gca().set_ylim((0.01, 0.1))
ylabel('k2')
xlabel('k1')
#plt.show()

figure(2)
x2, y2, prob2 = trace['k3']['k3'], trace['k3']['k4'],trace['k3']['prob']
x2, y2, prob2 = trace['k4']['k4'], trace['k4']['k3'],trace['k4']['prob']
scatter(x2, y2, s= 30, c= prob2)
#plt.gca().set_xlim((0, 0.1))
#plt.gca().set_ylim((0, 0.1))
xlabel('k3')
ylabel('k4')

1 Ответ

1 голос
/ 24 июня 2019

Для вашего примера, k3 и / или k4 в конечном итоге окажутся по крайней мере очень, очень близко к вашей верхней границе 1,0.Вы должны экономно использовать границы и иметь физическую значимость.То есть настоящая «проблема» заключается в том, что вы используете Parameters.add_many() и устанавливаете все параметры в пределах [0, 1].

В этом смысле подгонка сработала и даже сделала то, что вы просили;).Но если вы ослабите эти границы, вы, вероятно, получите лучший результат.

FWIW, вероятно, нет необходимости сначала решать с Nelder-Mead, а затем с leastsq (и если вы действительно хотите это сделать, вы должны правильно произнести Nelder0.

...