Exemplos de ajuste não-linear de parâmetros¶

Há vários pacotes de otimização disponíveis no julia
Aqui, vamos usar apenas o LsqFit e o Optim:
JuliaNLSolvers/LsqFit.jl: implementação do algoritmo de Levenberg-Maquardt especificamente para mínimos quadrados não-linear.
JuliaNLSolvers/Optim.jl: implementação de vários algoritmos de otimização, com qualquer função-objetivo, não necessariamente mínimos quadrados.

Outros pacotes¶

JuliaNLSolvers/LineSearches.jl: implementação de vários algoritmos do tipo busca linear (linesearch) (como o gradiente descendente, etc.).
jump-dev/JuMP.jl: para processos de modelagem mais complexos, com uma série de macros para construir um modelo para a otimização, e com interface para vários pacotes (em julia ou não) de otimização.
matthieugomez/LeastSquaresOptim.jl: também para a resolução de problemas de mínimos quadrados não lineares, via Levenberg-Marquardt ou Dogleg.
robertfeldt/BlackBoxOptim.jl: para otimizacões globais com métodos heurísticos/estocásticos, livres de derivada.
JuliaNLSolvers/NLsolve.jl: Para resolver sistemas de equações.
Veja mais em JuliaOpt

In [1]:

using Plots
using Random
using LsqFit
using Optim

Exemplo Michaelis-Menten de reação enzimática¶

Michaelis-Menten é um modelo para a taxa de reação enzimática (velocidade) $\nu$, dado por

$$ \nu = \frac{\nu_m t}{K_M + t} $$

Onde $t$ é a variável temporal, $\nu_m$ é uma taxa máxima de velodidade, e $K_M$ é um parâmetro de saturação da reação.
$\nu_m$ e $K_M$ são os parâmetros e, como de costume, denotamo-os por $\beta=(\nu_m, K_m)$.

In [2]:

function model(t, β)
    ν_m = β[1]
    K_M = β[2]
    v = (ν_m .* t) ./ (K_M .+ t)
    return v
end

ν_m = 0.3
K_M = 0.5
β = [ν_m K_M]
nothing

Obtendo uma amostra de dados¶

Perturbamos os dados exatos para simular imprecisões na coleta de dados

In [4]:

data_t = [0.03, 0.15, 0.3, 0.6, 1.3, 2.5, 3.5]
data_v = model(data_t, β) .+ 0.05*(rand(MersenneTwister(321), length(data_t)) .- 0.5)
plot(0:0.1:4, t -> model(t, β), label="taxa de reação", legend=:bottomright)
scatter!(data_t, data_v)
plot!(data_t, data_v, seriestype=:scatter, label="amostra com ruído")

Ajuste dos dados via LsqFit¶

Otimizando com o curve_fit do LsqFit
É necessário dar um "chute" inicial para os parâmetros.

In [5]:

β₀ = [0.5, 0.5]
fit = curve_fit(model, data_t, data_v, β₀)
β_fit = fit.param

2-element Vector{Float64}:
 0.27385430554237405
 0.3573547044880839

Visualizando o resultado¶

In [6]:

plot(0:0.1:4, t -> model(t, β), label="taxa de reação", legend=:bottomright)
plot!(data_t, data_v, seriestype=:scatter, label="amostra com ruído")
plot!(0:0.1:4, t -> model(t, β_fit), label="modelo ajustado", legend=:bottomright)

Ajuste dos dados via Optim¶

Aqui, usamos o método de gradiente descendente do Optim.
No Optim, precisamos explicitar a função custo.
Lembremos que queremos ajustar os parâmetros, ou seja, precisamos minimizar o custo em relação aos parâmetros.
Novamente, a função custo é da soma dos quadrados dos residuos, ou erro quadrático.

In [7]:

custo(β) = sum(abs2,data_v - model(data_t,β))

custo (generic function with 1 method)

In [9]:

resultado = optimize(custo, β₀, GradientDescent())

 * Status: success

 * Candidate solution
    Final objective value:     7.751687e-04

 * Found with
    Algorithm:     Gradient Descent

 * Convergence measures
    |x - x'|               = 7.31e-09 ≰ 0.0e+00
    |x - x'|/|x'|          = 2.05e-08 ≰ 0.0e+00
    |f(x) - f(x')|         = 1.55e-16 ≰ 0.0e+00
    |f(x) - f(x')|/|f(x')| = 2.01e-13 ≰ 0.0e+00
    |g(x)|                 = 7.10e-09 ≤ 1.0e-08

 * Work counters
    Seconds run:   0  (vs limit Inf)
    Iterations:    72
    f(x) calls:    220
    ∇f(x) calls:   220

In [10]:

β_optim = Optim.minimizer(resultado)

2-element Vector{Float64}:
 0.2738543232663811
 0.3573548106304591

Visualizando o resultado¶

In [11]:

plot(0:0.1:4, t -> model(t, β), label="taxa de reação", legend=:bottomright)
plot!(data_t, data_v, seriestype=:scatter, label="amostra com ruído")
plot!(0:0.1:4, t -> model(t, β_fit), label="modelo ajustado via Levenberg-Marquardt do LsqFit", legend=:bottomright)
plot!(0:0.1:4, t -> model(t, β_optim), seriestype=:scatter, markershape=:xcross, markersize=2,
        label="modelo ajustado via gradiente descendente do Optim", legend=:bottomright)

Outros métodos de otimização¶

No Optim, podemos usar diversos outros métodos, basta substituir o argumento GradientDescent() pelo método desejado.
Podemos usar métodos que usam o gradiente: ConjugateGradient(), GradientDescent(), BFGS() e LBFSG().
Métodos que usam a Hessiana: Newton(), NewtonTrustRegion(), IPNewton() - interior point.
Métodos livres de gradiente: NelderMead(), SimulatedAnnealing(), SAMIN() - simulated annealing with bounds, ParticleSwarm().
Cada construtor acima aceita diversos parâmetros.
Por default, diferenças finitas são usadas para calcular o gradiente e a Hessiana, quando necessário. Mas pode se passar a derivada expliciamente (passando os argumentos opcionais f! e g!) ou se usar derivação automática (passando-se o argumento autodiff, e.g. autodiff = :forward).

In [24]:

optimize(custo, β₀, LBFGS())

 * Status: success

 * Candidate solution
    Final objective value:     7.751687e-04

 * Found with
    Algorithm:     L-BFGS

 * Convergence measures
    |x - x'|               = 6.46e-08 ≰ 0.0e+00
    |x - x'|/|x'|          = 1.81e-07 ≰ 0.0e+00
    |f(x) - f(x')|         = 2.28e-16 ≰ 0.0e+00
    |f(x) - f(x')|/|f(x')| = 2.94e-13 ≰ 0.0e+00
    |g(x)|                 = 3.31e-12 ≤ 1.0e-08

 * Work counters
    Seconds run:   1  (vs limit Inf)
    Iterations:    9
    f(x) calls:    25
    ∇f(x) calls:   25

In [25]:

optimize(custo, β₀, Newton(), autodiff = :forward)

 * Status: success

 * Candidate solution
    Final objective value:     7.751687e-04

 * Found with
    Algorithm:     Newton's Method

 * Convergence measures
    |x - x'|               = 2.34e-06 ≰ 0.0e+00
    |x - x'|/|x'|          = 6.56e-06 ≰ 0.0e+00
    |f(x) - f(x')|         = 1.11e-11 ≰ 0.0e+00
    |f(x) - f(x')|/|f(x')| = 1.43e-08 ≰ 0.0e+00
    |g(x)|                 = 2.39e-11 ≤ 1.0e-08

 * Work counters
    Seconds run:   1  (vs limit Inf)
    Iterations:    6
    f(x) calls:    20
    ∇f(x) calls:   20
    ∇²f(x) calls:  6

Dados espúrios (outliers)¶

É comum, por dificuldades no processo de amostragem ou de tratamento dos dados, termos dados "espúrios", i.e. claramente fora do padrão.
São, também, denominados de outliers.
Podemos tentar filtrar os dados espúrios da amostra usando técnicas de estatística.
E/Ou podemos reduzir o impacto dos dados espúrios com funções de custo ligeiramente modificadas.

Dados contaminados¶

Vamos modificar a amostra anterior substituindo um dos dados por um dado espúrio.

In [12]:

data_v_esp = copy(data_v)
data_v_esp[5] += 0.15
nothing

In [13]:

plot(0:0.1:4, t -> model(t, β), label="taxa de reação", legend=:bottomright)
plot!(data_t, data_v_esp, seriestype=:scatter, label="amostra com ruído espúrio")

Ajuste¶

In [14]:

custo_esp(β) = sum(abs2,data_v_esp - model(data_t,β))
resultado_esp = optimize(custo_esp, β₀, GradientDescent())
β_optim_esp = Optim.minimizer(resultado_esp)

2-element Vector{Float64}:
 0.3150999020168445
 0.3607669696333052

In [15]:

plot(0:0.1:4, t -> model(t, β), label="taxa de reação", legend=:bottomright)
plot!(data_t, data_v_esp, seriestype=:scatter, label="amostra com ruído espúrio")
plot!(0:0.1:4, t -> model(t, β_optim_esp), linestyle=:dash,
        label="modelo ajustado via gradiente descendente do Optim", legend=:bottomright)

Funções de custo modificadas¶

A ideia, aqui, é modificar a função quadrática de custo para uma em que resíduos relativamente grandes causem menos influência.
Vários tipos de modificações costumam ser usadas.

In [16]:

mods_custo = Dict(
    :linear => (r -> @. r),
    :soft_l1 => (r -> @. 2*((1 + abs(r))^(1/2) - 1)),
    :huber => (r -> @. min(r, 2*((1 + abs(r))^(1/2) - 1))),
    :cauchy => (r -> @. log.(1 .+ abs(r))),
    :arctan => (r -> @. atan(r))
)

Dict{Symbol, Function} with 5 entries:
  :cauchy  => #32
  :arctan  => #33
  :soft_l1 => #30
  :huber   => #31
  :linear  => #29

In [17]:

plot(title="Diferentes modificadores para a função de custo", titlefont=10)
for (k,func) in mods_custo
    plot!(0:0.1:1.5, func, label="função modificadora $k", 
        legendfont=6, legend=:topleft, subplot=1)
end
plot!()

In [18]:

sq(e) = sum(abs2, e)

sq (generic function with 1 method)

In [19]:

plot(title="Diferentes modificadores para a função de custo", titlefont=10)
for (k,func) in mods_custo
    plot!(0:0.1:1.5, sq∘func, label="função modificadora $k", 
        legendfont=6, legend=:topleft, subplot=1)
end
plot!()

Ampliando o efeito dos modificadores¶

O efeito dos diferents modificadores pode ser amplificado com uma função de espalhamento:

$$ s \mapsto \lambda s $$

In [20]:

λ(s) = 20*s

λ (generic function with 1 method)

In [21]:

plot(title="Diferentes modificadores para a função de custo", titlefont=10)
for (k,func) in mods_custo
    plot!(0:0.1:1.5, sq∘λ∘func, label="função modificadora $k", 
        legendfont=6, legend=:topleft, subplot=1)
end
plot!()

Função de custo modificada¶

Assim, a função de custo é a composição

$$ \text{função quadrática} \circ \text{modificador} \circ \text{espalhamento} \circ \text{resíduos}. $$

Fazendo os ajustes¶

In [22]:

residuos_esp(β) = data_v_esp - model(data_t,β)
resultados_esp = 
    Dict(k => optimize(sq∘func∘λ∘residuos_esp, β₀, GradientDescent())
        for (k,func) in mods_custo
        )
nothing

Visualizando os resultados¶

In [23]:

plot(0:0.1:4, t -> model(t, β), label="taxa de reação", legend=:bottomright)
plot!(data_t, data_v_esp, seriestype=:scatter, label="amostra com ruído espúrio")
for k in keys(mods_custo)
    plot!(0:0.1:4, t -> model(t, Optim.minimizer(resultados_esp[k])), linestyle=:dash,
        label="modelo ajustado com modificador de custo $k", legend=:bottomright)
end
plot!(title="Ajustes com diferentes modificadores da função de custo", titlefont=10)

Exercícios¶

Refaça os exemplos do caderno sobre Modelos redutíveis ao caso linear nos parâmetros e aplicações como problemas de otimização não-linear, sem usar as mudanças de variáveis feitas lá.