Оценка параметров множественных наборов данных в юлиях

Question

Я искал и не смог найти прямой способ использования оценки параметра diffrentialEquations в julia для подбора нескольких наборов данных.Итак, допустим, у нас есть это простое дифференциальное уравнение с двумя параметрами:

f1 = function (du,u,p,t)
  du[1] = - p[1]*p[2] * u[1]
end

У нас есть экспериментальные наборы данных u [1] vs t.Каждый набор данных имеет различное значение p [2] и / или разные начальные условия.p [1] - параметр, который мы хотим оценить.Я могу сделать это, решив дифференциальное уравнение в цикле for, который перебирает различные начальные условия и значения p [2], сохраняя решения в массиве и создавая функцию потерь по экспериментальным данным.Интересно, есть ли способ сделать это в меньшем количестве строк кода, используя, например, DiffEqBase.problem_new_parameters для установки условий каждого набора данных.Это очень распространенная ситуация при подгонке моделей к экспериментальным данным, но я не смог найти хороший пример в документации.

Заранее спасибо,

С уважением

РЕДАКТИРОВАТЬ 1

Вышеуказанный случай является лишь упрощенным примером.Чтобы сделать это на практике, мы могли бы создать некоторые поддельные экспериментальные данные из следующего кода:

using DifferentialEquations

# ODE function
f1 = function (du,u,p,t)
  du[1] = - p[1]*p[2] * u[1]
end

# Initial conditions and parameter values.
# p1 is the parameter to be estimated.
# p2 and u0 are experimental parameters known for each dataset.
u0 = [1.,2.]
p1 = .3
p2 = [.1,.2]
tspan = (0.,10.)
t=linspace(0,10,5)


# Creating data for 1st experimental dataset.
## Experimental conditions: u0 = 1. ; p[2] = .1
prob1 = ODEProblem(f1,[u0[1]],tspan,[p1,p2[1]])
sol1=solve(prob1,Tsit5(),saveat=t)

# Creating data for 2nd experimental dataset.
## Experimental conditions: u0 = 1. ; p[2] = .2
prob2 = ODEProblem(f1,[u0[1]],tspan,[p1,p2[2]])
sol2=solve(prob2,Tsit5(),saveat=t)

# Creating data for 3rd experimental dataset.
## Experimental conditions: u0 = 2. ; p[2] = .1
prob3 = ODEProblem(f1,[u0[2]],tspan,[p1,p2[1]])
sol3=solve(prob3,Tsit5(),saveat=t)

sol1, sol2 и sol3 теперь являются нашими экспериментальными данными, каждый набор данных использует различную комбинацию начальных условий и p [2] (которая представляет некоторую экспериментальную переменную (например, температуру, расход ...)

Цель состоит в том, чтобы оценить значение p [1], используя экспериментальные данные sol1, sol2 и sol3, позволяющие DiffEqBase.problem_new_parameters илидругая альтернатива итерации по экспериментальным условиям.

Answer 1

Что вы можете сделать, это создать проблему MonteCarloProblem, которая решает сразу все три проблемы:

function prob_func(prob,i,repeat)
  i < 3 ? u0 = [1.0] : u0 = [2.0]
  i == 2 ? p = (prob.p[1],0.2) : p = (prob.p[1],0.1)
  ODEProblem{true}(f1,u0,(0.0,10.0),p)
end
prob = MonteCarloProblem(prob1,prob_func = prob_func)

@time sol = solve(prob,Tsit5(),saveat=t,parallel_type = :none,
                  num_monte = 3)

Затем создайте функцию потерь, которая сравнивает каждое из решений с 3 наборами данных и складывает их потери вместе..

loss1 = L2Loss(t,data1)
loss2 = L2Loss(t,data2)
loss3 = L2Loss(t,data3)
loss(sol) = loss1(sol[1]) + loss2(sol[2]) + loss3(sol[3])

Наконец, вам необходимо рассказать, как соотнести параметры оптимизации с проблемой, которую они решают.Здесь наша проблема MonteCarloProblem содержит prob, которую она вытягивает p[1] из любого случая, когда она создает проблему.Значение, которое мы хотим оптимизировать, - это p[1], поэтому:

function my_problem_new_parameters(prob,p)
  prob.prob.p[1] = p[1]
  prob
end

Теперь наша цель состоит именно из этих частей:

obj = build_loss_objective(prob,Tsit5(),loss,
                           prob_generator = my_problem_new_parameters,
                           num_monte = 3,
                           saveat = t)

Теперь давайте бросим это в Optim.jl'sМетод Брента:

using Optim
res = optimize(obj,0.0,1.0)

Results of Optimization Algorithm
 * Algorithm: Brent's Method
 * Search Interval: [0.000000, 1.000000]
 * Minimizer: 3.000000e-01
 * Minimum: 2.004680e-20
 * Iterations: 10
 * Convergence: max(|x - x_upper|, |x - x_lower|) <= 2*(1.5e-08*|x|+2.2e-16): true
 * Objective Function Calls: 11

Обнаружено, что наилучшим общим значением является 0.3, который является параметром, который мы использовали для генерации данных.

Вот полный код:

using DifferentialEquations

# ODE function
f1 = function (du,u,p,t)
  du[1] = - p[1]*p[2] * u[1]
end

# Initial conditions and parameter values.
# p1 is the parameter to be estimated.
# p2 and u0 are experimental parameters known for each dataset.
u0 = [1.,2.]
p1 = .3
p2 = [.1,.2]
tspan = (0.,10.)
t=linspace(0,10,5)

# Creating data for 1st experimental dataset.
## Experimental conditions: u0 = 1. ; p[2] = .1
prob1 = ODEProblem(f1,[u0[1]],tspan,[p1,p2[1]])
sol1=solve(prob1,Tsit5(),saveat=t)
data1 = Array(sol1)

# Creating data for 2nd experimental dataset.
## Experimental conditions: u0 = 1. ; p[2] = .2
prob2 = ODEProblem(f1,[u0[1]],tspan,[p1,p2[2]])
sol2=solve(prob2,Tsit5(),saveat=t)
data2 = Array(sol2)

# Creating data for 3rd experimental dataset.
## Experimental conditions: u0 = 2. ; p[2] = .1
prob3 = ODEProblem(f1,[u0[2]],tspan,[p1,p2[1]])
sol3=solve(prob3,Tsit5(),saveat=t)
data3 = Array(sol3)

function prob_func(prob,i,repeat)
  i < 3 ? u0 = [1.0] : u0 = [2.0]
  i == 2 ? p = (prob.p[1],0.2) : p = (prob.p[1],0.1)
  ODEProblem{true}(f1,u0,(0.0,10.0),p)
end
prob = MonteCarloProblem(prob1,prob_func = prob_func)

# Just to show what this looks like
sol = solve(prob,Tsit5(),saveat=t,parallel_type = :none,
            num_monte = 3)

loss1 = L2Loss(t,data1)
loss2 = L2Loss(t,data2)
loss3 = L2Loss(t,data3)
loss(sol) = loss1(sol[1]) + loss2(sol[2]) + loss3(sol[3])

function my_problem_new_parameters(prob,p)
  prob.prob.p[1] = p[1]
  prob
end
obj = build_loss_objective(prob,Tsit5(),loss,
                           prob_generator = my_problem_new_parameters,
                           num_monte = 3,
                           saveat = t)

using Optim
res = optimize(obj,0.0,1.0)