2

非線形最適化に lmfit python パッケージを使用しています (url: http://lmfit.github.io/lmfit-py/ )。最小二乗法を使用する場合、ヤコビ関数を渡すことが可能かどうかを知りたいですか? はいの場合、最小限の例を提供していただけますか?

ありがとうございました！カブラルプ

PS: コードは次のとおりです。

# f(t,g,p) = dg_dt(t,g,p) = R*(c^h/(c^h+K^h))-l*g
# returns rhs of an ODE (dg_dt)
def hill_1g1c(t, g_in, p_in):
    R = p_in['R'].value
    K = p_in['K'].value
    h = p_in['h'].value
    l = p_in['l'].value

    dg_dt = R*((c_int(t)**h)/((c_int(t)**h)+(K**h))) - l*g_in
    return dg_dt

# f_deriv(t,g,p)
# is intended to return the derivatives of f with respect to 4 parameter
def hill_1g1c_jac(p_in, y_in, dt, t_max, g_data, g_err):
    t=1
    R = p_in['R'].value
    K = p_in['K'].value
    h = p_in['h'].value
    l = p_in['l'].value
    dg_dR = (c_int(t)**h) / (c_int(t)**h + K**h) - l * 1
    dg_dK = (-1 * R * c_int(t)**h * h * k**(h-1)) / ((c_int(t)**h + K**h)**2) - l * 1
    dg_dh = (-1 * R * c_int(t)**h * k**h * (np.log(k) - np.log(c_int(t)))) / ((c_int(t)**h + K**h)**2) - l * 1
    dg_dl = -y_in - l * 1

    return np.array([dg_dR, dg_dK, dg_dh, dg_dl])

# y = ODE_solve(y0,p,dt, t_max) >>>> wrapper around ode.integrate, returns array of g
def ODE_solve(y0, p_in, dt, t_max):
    t = [0]
    y = [y0]
    r.set_initial_value(y0, t=0.0)
    r.set_f_params(p_in)
    while r.successful() and r.t < t_max:
         r.integrate(r.t+dt)
         t.append(r.t)
         y.append(r.y)
    return np.array(y)

# weighted least squares, objective function to be minimised
def ODE_wres(p_in, y0, dt, t_max, g_data, g_err):
    g_extended = ODE_solve(y0, p_in, dt, t_max)
    g_model = g_extended[-25:]/g_extended[-25]
    weighted_residuals = (g_data - g_model)/(g_err + 0.00000001)
    return weighted_residuals


# specs for inegrate.ode
y0 = 1
t0 = 0
r = integrate.ode(hill_1g1c).set_integrator('vode', with_jacobian=False)
t_stim = 15
t_max = t_stim + 24
t_plus = 5
dt = 1
t_extended = np.linspace(0,t_max+t_plus,t_max+t_plus+1)

# set history of all inputs to 1
c_history = [1 for val in range(t_stim)]

# data (is read in from text file)
g_data = Y_data[:,i]
# error in g
g_err = Y_error[:,i]
# input c
c_data = Y_data[:,k]
# interpolation of c, contains history (=1) and future (=endval)
c_future = [c_data[-1] for val in range(t_plus)]
c_extended = np.hstack((c_history, c_data, c_future))
c_int = interp1d(t_extended, c_extended, kind='linear')

# initial parameter vector
R_ini = random.uniform(0.01, 500.0)
K_ini = random.uniform(0.01, 20.0)
h_ini = random.uniform(-4.0, 4.0)
l_ini = random.uniform(0.07, 7.0)

p_ini = Parameters()
p_ini.add('R', value= R_ini, min=0.01, max=500)
p_ini.add('K', value= K_ini, min=0.01, max=20)
p_ini.add('h', value= h_ini, min=-4, max=4)
p_ini.add('l', value= l_ini, min=0.07, max=7.0)

res_ini = ODE_wres(p_ini, y0, dt, t_max, g_data, g_err)
chisqr_ini = np.sum(res_ini**2)

# optimise
lmsol = Minimizer(ODE_wres, p_ini, fcn_args=(y0, dt, t_max, g_data, g_err))
lmsol.leastsq(Dfun=hill_1g1c_jac, col_deriv=True)

PPS: github でこの貴重な例を見つけました: https://github.com/lmfit/lmfit-py/blob/master/examples/example_derivfunc.py

注意: ヤコビ関数を lmft.leastsq に渡すことができた後、私のテスト ケースでは、lmfit によって返された最適化された解が真の解に収束していないことに気付きました。ただし、実際の scipy.optimize.leastq 関数 (lmfit によって呼び出される) を使用すると、すべて正常に機能しました。つまり、返された解は、ヤコビアンも含めて収束しました。ヤコビ関数を指定すると lmfit.leastsq が適切に機能しないと言っているわけではありませんが、この場合は注意して扱うことをお勧めします。これまでのところ、この原因を詳しく調べる時間はありませんでした。

4

1 に答える 1

3

Dfunメソッドへの引数としてヤコビアンを計算する関数を渡すことができますMinimizer.leastsq: http://lmfit.github.io/lmfit-py/fitting.html#leastsq

デフォルトでは、 for に渡された関数は、Dfunパラメーターと同じ行数の配列を返し、各行は適合される各パラメーターに関する導関数を返します。パラメータが正しい順序で処理されるように、 Parametersオブジェクトを使用してパラメータを指定してください。関係ないかもしれませんが、これは必要な IIRC だと思います。

于 2014-01-02T22:08:21.467 に答える