Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Scalarize input #202

Open
wants to merge 5 commits into
base: master
Choose a base branch
from
Open

Scalarize input #202

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
046b6d7
Relax tests
gustaphe Oct 23, 2021
e09353a
Make model accept scalar arguments.
gustaphe Oct 24, 2021
00ec20d
Fix 'scalar' multivariate fit
gustaphe Oct 25, 2021
d0339fb
Update README
gustaphe Oct 25, 2021
3cfaa7f
Update docs
gustaphe Oct 25, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
51 changes: 31 additions & 20 deletions README.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,21 +14,19 @@ There are top-level methods `curve_fit()` and `margin_error()` that are useful f
using LsqFit

# a two-parameter exponential model
# x: array of independent variables
# x: independent variable
# p: array of model parameters
# model(x, p) will accept the full data set as the first argument `x`.
# This means that we need to write our model function so it applies
# the model to the full dataset. We use `@.` to apply the calculations
# across all rows.
@. model(x, p) = p[1]*exp(-x*p[2])
# model(x, p) will accept the independent variable as the first argument `x`.
model(x, p) = p[1]*exp(-x*p[2])
```
The function applies the per observation function `p[1]*exp(-x[i]*p[2])` to the full dataset in `x`, with `i` denoting an observation row. We simulate some data and chose our "true" parameters.
We simulate some data and chose our "true" parameters.

```julia
# some example data
# xdata: independent variables
# ydata: dependent variable
xdata = range(0, stop=10, length=20)
ydata = model(xdata, [1.0 2.0]) + 0.01*randn(length(xdata))
ydata = model.(xdata, Ref([1.0 2.0])) + 0.01*randn(length(xdata))
p0 = [0.5, 0.5]
```
Now, we're ready to fit the model.
Expand Down Expand Up @@ -56,21 +54,31 @@ confidence_inter = confidence_interval(fit, 0.05)

# The finite difference method is used above to approximate the Jacobian.
# Alternatively, a function which calculates it exactly can be supplied instead.
# It should return a vector `J_i = df/d(p_i)`
function jacobian_model(x,p)
J = Array{Float64}(undef, length(x), length(p))
@. J[:,1] = exp(-x*p[2]) #dmodel/dp[1]
@. @views J[:,2] = -x*p[1]*J[:,1] #dmodel/dp[2], thanks to @views we don't allocate memory for the J[:,1] slice
J = Array{Float64}(undef, length(p))
J[1] = exp(-x*p[2]) #dmodel/dp[1]
J[2] = -x*p[1]*J[1] #dmodel/dp[2]
J
end
fit = curve_fit(model, jacobian_model, xdata, ydata, p0)
```

Multivariate regression
-----------------------
There's nothing inherently different if there are more than one variable entering the problem. We just need to specify the columns appropriately in our model specification:
There's nothing inherently different if there are more than one variable
entering the problem. The argument `x` must just be a vector:
```julia
multimodel(x, p) = p[1]*exp(-x[1]*p[2]+x[2]*p[3])
```
and the `x` argument to the multivariate `curve_fit` function should be a
vector of such input vectors.
```julia
@. multimodel(x, p) = p[1]*exp(-x[:, 1]*p[2]+x[:, 2]*p[3])
x = [randn(2) for _ in 1:50]
y = multimodel.(x, Ref([1.0, 2.0, 1.0])) + randn(50)
fit = curve_fit(multimodel, x, y, [1.0, 1.0, 1.0])
```

Evaluating the Jacobian and using automatic differentiation
-------------------------
The default is to calculate the Jacobian using a central finite differences scheme if no Jacobian function is provided. The default is to use central differences because it can be more accurate than forward finite differences, but at the expense of computational cost. It is possible to switch to forward finite differences, like MINPACK uses for example, by specifying `autodiff=:finiteforward`:
Expand All @@ -85,15 +93,18 @@ Here, you have to be careful not to manually restrict any types in your code to,

In-place model and Jacobian
-------------------------
It is possible to either use an in-place model, or an in-place model *and* an in-place Jacobian. It might be pertinent to use this feature when `curve_fit` is slow, or consumes a lot of memory.
It is possible to either use an in-place model, or an in-place model *and* an
in-place Jacobian. It might be pertinent to use this feature when `curve_fit`
is slow, or consumes a lot of memory. An in-place model must operate on the
entire input vector, not individual elements.
```
model_inplace(F, x, p) = (@. F = p[1] * exp(-x * p[2]))
model_inplace!(F, x, p) = (@. F = p[1] * exp(-x * p[2]))

function jacobian_inplace(J::Array{Float64,2},x,p)
function jacobian_inplace!(J::Array{Float64,2},x,p)
@. J[:,1] = exp(-x*p[2])
@. @views J[:,2] = -x*p[1]*J[:,1]
end
fit = curve_fit(model_inplace, jacobian_inplace, xdata, ydata, p0; inplace = true)
fit = curve_fit(model_inplace!, jacobian_inplace!, xdata, ydata, p0; inplace = true)
```

Geodesic acceleration
Expand Down Expand Up @@ -161,7 +172,7 @@ This performs a fit using a non-linear iteration to minimize the (weighted) resi

`sigma = stderror(fit; atol, rtol)`:

* `fit`: result of curve_fit (a `LsqFitResult` type)
* `fit`: result of `curve_fit` (a `LsqFitResult` type)
* `atol`: absolute tolerance for negativity check
* `rtol`: relative tolerance for negativity check

Expand All @@ -171,7 +182,7 @@ If no weights are provided for the fits, the variance is estimated from the mean

`margin_of_error = margin_error(fit, alpha=0.05; atol, rtol)`:

* `fit`: result of curve_fit (a `LsqFitResult` type)
* `fit`: result of `curve_fit` (a `LsqFitResult` type)
* `alpha`: significance level
* `atol`: absolute tolerance for negativity check
* `rtol`: relative tolerance for negativity check
Expand All @@ -180,7 +191,7 @@ This returns the product of standard error and critical value of each parameter

`confidence_interval = confidence_interval(fit, alpha=0.05; atol, rtol)`:

* `fit`: result of curve_fit (a `LsqFitResult` type)
* `fit`: result of `curve_fit` (a `LsqFitResult` type)
* `alpha`: significance level
* `atol`: absolute tolerance for negativity check
* `rtol`: relative tolerance for negativity check
Expand Down
9 changes: 5 additions & 4 deletions docs/src/getting_started.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -16,9 +16,9 @@ m(t, \boldsymbol{p}) = p_1 \exp(-p_2 t)
```

```julia
julia> # t: array of independent variable
julia> # t: independent variable
julia> # p: array of model parameters
julia> model(t, p) = p[1] * exp.(-p[2] * t)
julia> model(t, p) = p[1] * exp(-p[2] * t)
```

For illustration purpose, we generate some fake data.
Expand All @@ -27,7 +27,7 @@ For illustration purpose, we generate some fake data.
julia> # tdata: data of independent variable
julia> # ydata: data of dependent variable
julia> tdata = range(0, stop=10, length=20)
julia> ydata = model(tdata, [1.0 2.0]) + 0.01*randn(length(tdata))
julia> ydata = model.(tdata, Ref([1.0 2.0])) + 0.01*randn(length(tdata))
```

Before fitting the data, we also need a initial value of parameters for `curve_fit()`.
Expand All @@ -46,7 +46,8 @@ julia> param = fit.param
2.0735
```

`LsqFit.jl` also provides functions to examinep0 = [0.5, 0.5] the goodness of fit. `estimate_covar(fit)` computes the estimated covariance matrix.
`LsqFit.jl` also provides functions to examine the goodness of fit.
`estimate_covar(fit)` computes the estimated covariance matrix.

```Julia
julia> cov = estimate_covar(fit)
Expand Down
12 changes: 6 additions & 6 deletions docs/src/tutorial.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -31,9 +31,9 @@ Y_i = \gamma_1 \exp(\gamma_2 t_i) + \epsilon_i
To fit data using `LsqFit.jl`, pass the defined model function (`m`), data (`tdata` and `ydata`) and the initial parameter value (`p0`) to `curve_fit()`. For now, `LsqFit.jl` only supports the Levenberg Marquardt algorithm.

```julia
julia> # t: array of independent variables
julia> # t: independent variable
julia> # p: array of model parameters
julia> m(t, p) = p[1] * exp.(p[2] * t)
julia> m(t, p) = p[1] * exp(p[2] * t)
julia> p0 = [0.5, 0.5]
julia> fit = curve_fit(m, tdata, ydata, p0)
```
Expand Down Expand Up @@ -80,9 +80,9 @@ By default, the finite differences is used (see [NLSolversBase.jl](https://githu

```Julia
function j_m(t,p)
J = Array{Float64}(length(t),length(p))
J[:,1] = exp.(p[2] .* t) #df/dp[1]
J[:,2] = t .* p[1] .* J[:,1] #df/dp[2]
J = Array{Float64}(length(p))
J[1] = exp(p[2] * t) #df/dp[1]
J[2] = t * p[1] * J[1] #df/dp[2]
J
end

Expand Down Expand Up @@ -267,7 +267,7 @@ The algorithm in `LsqFit.jl` will then provide a least squares solution $\boldsy
In `LsqFit.jl`, the residual function passed to `levenberg_marquardt()` is in different format, if the weight is a vector:

```julia
r(p) = sqrt.(wt) .* ( model(xpts, p) - ydata )
r(p) = sqrt.(wt) .* ( model.(xpts, Ref(p)) - ydata )
lmfit(r, g, p0, wt; kwargs...)
```

Expand Down
33 changes: 18 additions & 15 deletions src/curve_fit.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -16,14 +16,17 @@ StatsBase.residuals(lfr::LsqFitResult) = lfr.resid
mse(lfr::LsqFitResult) = rss(lfr)/dof(lfr)

function check_data_health(xdata, ydata)
if any(ismissing, xdata) || any(ismissing, ydata)
error("Data contains `missing` values and a fit cannot be performed")
if is_unhealthy(xdata)
error("x data contains `missing`, `Inf` or `NaN` values and a fit cannot be performed")
end
if any(isinf, xdata) || any(isinf, ydata) || any(isnan, xdata) || any(isnan, ydata)
error("Data contains `Inf` or `NaN` values and a fit cannot be performed")
if is_unhealthy(ydata)
error("y data contains `missing`, `Inf` or `NaN` values and a fit cannot be performed")
end
end

is_unhealthy(x) = ismissing(x) || isinf(x) || isnan(x)
is_unhealthy(x::AbstractArray) = any(is_unhealthy, x)

# provide a method for those who have their own Jacobian function
function lmfit(f, g, p0::AbstractArray, wt::AbstractArray; kwargs...)
r = f(p0)
Expand All @@ -47,7 +50,7 @@ function lmfit(f, p0::AbstractArray, wt::AbstractArray; autodiff = :finite, kwar
# this is a convenience function for the curve_fit() methods
# which assume f(p) is the cost functionj i.e. the residual of a
# model where
# model(xpts, params...) = ydata + error (noise)
# model(x, params...) = y + error (noise)

# this minimizes f(p) using a least squares sum of squared error:
# rss = sum(f(p)^2)
Expand Down Expand Up @@ -97,7 +100,7 @@ model(x, p) = p[1]*exp.(-x.*p[2])
# xdata: independent variables
# ydata: dependent variable
xdata = range(0, stop=10, length=20)
ydata = model(xdata, [1.0 2.0]) + 0.01*randn(length(xdata))
ydata = model.(xdata, Ref([1.0 2.0])) .+ 0.01*randn(length(xdata))
p0 = [0.5, 0.5]

fit = curve_fit(model, xdata, ydata, p0)
Expand All @@ -114,7 +117,7 @@ function curve_fit(model, xdata::AbstractArray, ydata::AbstractArray, p0::Abstra
f! = (F,p) -> (model(F,xdata,p); @. F = F - ydata)
lmfit(f!, p0, T[], ydata; kwargs...)
else
f = (p) -> model(xdata, p) - ydata
f = (p) -> model.(xdata, Ref(p)) - ydata
lmfit(f, p0, T[]; kwargs...)
end
end
Expand All @@ -130,8 +133,8 @@ function curve_fit(model, jacobian_model,
g! = (G,p) -> jacobian_model(G, xdata, p)
lmfit(f!, g!, p0, T[], copy(ydata); kwargs...)
else
f = (p) -> model(xdata, p) - ydata
g = (p) -> jacobian_model(xdata, p)
f = (p) -> model.(xdata, Ref(p)) - ydata
g = (p) -> vcat(reshape.(jacobian_model.(xdata, Ref(p)), 1, :)...)
lmfit(f, g, p0, T[]; kwargs...)
end
end
Expand All @@ -146,7 +149,7 @@ function curve_fit(model, xdata::AbstractArray, ydata::AbstractArray, wt::Abstra
f! = (F,p) -> (model(F,xdata,p); @. F = u*(F - ydata))
lmfit(f!, p0, wt, ydata; kwargs...)
else
f = (p) -> u .* ( model(xdata, p) - ydata )
f = (p) -> u .* ( model.(xdata, Ref(p)) - ydata )
lmfit(f,p0,wt; kwargs...)
end
end
Expand All @@ -161,8 +164,8 @@ function curve_fit(model, jacobian_model,
g! = (G,p) -> (jacobian_model(G, xdata, p); @. G = u*G )
lmfit(f!, g!, p0, wt, ydata; kwargs...)
else
f = (p) -> u .* ( model(xdata, p) - ydata )
g = (p) -> u .* ( jacobian_model(xdata, p) )
f = (p) -> u .* ( model.(xdata, Ref(p)) - ydata )
g = (p) -> u .* ( vcat(reshape.(jacobian_model.(xdata, Ref(p)), 1, :)...) )
lmfit(f, g, p0, wt; kwargs...)
end
end
Expand All @@ -179,7 +182,7 @@ function curve_fit(model, xdata::AbstractArray, ydata::AbstractArray, wt::Abstra
# This requires the matrix to be positive definite
u = cholesky(wt).U

f(p) = u * ( model(xdata, p) - ydata )
f(p) = u * ( model.(xdata, Ref(p)) - ydata )
lmfit(f,p0,wt; kwargs...)
end

Expand All @@ -189,8 +192,8 @@ function curve_fit(model, jacobian_model,

u = cholesky(wt).U

f(p) = u * ( model(xdata, p) - ydata )
g(p) = u * ( jacobian_model(xdata, p) )
f(p) = u * ( model.(xdata, Ref(p)) - ydata )
g(p) = u * ( vcat(reshape.(jacobian_model.(xdata, Ref(p)), 1, :)...) )
lmfit(f, g, p0, wt; kwargs...)
end

Expand Down
37 changes: 18 additions & 19 deletions test/curve_fit.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,63 +1,62 @@
let
@testset "curve_fit" begin
# before testing the model, check whether missing/null data is rejected
tdata = [rand(1:10, 5)..., missing]
@test_throws ErrorException("Data contains `missing` values and a fit cannot be performed") LsqFit.check_data_health(tdata, tdata)
@test_throws ErrorException("x data contains `missing`, `Inf` or `NaN` values and a fit cannot be performed") LsqFit.check_data_health(tdata, tdata)
tdata = [rand(1:10, 5)..., Inf]
@test_throws ErrorException("Data contains `Inf` or `NaN` values and a fit cannot be performed") LsqFit.check_data_health(tdata, tdata)
@test_throws ErrorException("x data contains `missing`, `Inf` or `NaN` values and a fit cannot be performed") LsqFit.check_data_health(tdata, tdata)
tdata = [rand(1:10, 5)..., NaN]
@test_throws ErrorException("Data contains `Inf` or `NaN` values and a fit cannot be performed") LsqFit.check_data_health(tdata, tdata)
@test_throws ErrorException("x data contains `missing`, `Inf` or `NaN` values and a fit cannot be performed") LsqFit.check_data_health(tdata, tdata)

# fitting noisy data to an exponential model
# TODO: Change to `.-x` when 0.5 support is dropped
model(x, p) = p[1] .* exp.(-x .* p[2])
model(x, p) = p[1] * exp(-x * p[2])

# some example data
Random.seed!(12345)
xdata = range(0, stop=10, length=20)
ydata = model(xdata, [1.0, 2.0]) + 0.01*randn(length(xdata))
ydata = model.(xdata, Ref([1.0, 2.0])) + 0.01*randn(length(xdata))
p0 = [0.5, 0.5]

for ad in (:finite, :forward, :forwarddiff)
fit = curve_fit(model, xdata, ydata, p0; autodiff = ad)
@assert norm(fit.param - [1.0, 2.0]) < 0.05
@test norm(fit.param - [1.0, 2.0]) < 0.05
@test fit.converged

# can also get error estimates on the fit parameters
errors = margin_error(fit, 0.1)
@assert norm(errors - [0.017, 0.075]) < 0.01
@test norm(errors - [0.017, 0.075]) < 0.015
end
# if your model is differentiable, it can be faster and/or more accurate
# to supply your own jacobian instead of using the finite difference
function jacobian_model(x,p)
J = Array{Float64}(undef, length(x), length(p))
J[:,1] = exp.(-x.*p[2]) #dmodel/dp[1]
J[:,2] = -x.*p[1].*J[:,1] #dmodel/dp[2]
J = Array{Float64}(undef, length(p))
J[1] = exp(-x*p[2]) #dmodel/dp[1]
J[2] = -x*p[1]*J[1] #dmodel/dp[2]
J
end
jacobian_fit = curve_fit(model, jacobian_model, xdata, ydata, p0)
@assert norm(jacobian_fit.param - [1.0, 2.0]) < 0.05
@test norm(jacobian_fit.param - [1.0, 2.0]) < 0.05
@test jacobian_fit.converged

# some example data
yvars = 1e-6*rand(length(xdata))
ydata = model(xdata, [1.0, 2.0]) + sqrt.(yvars) .* randn(length(xdata))
ydata = model.(xdata, Ref([1.0, 2.0])) + sqrt.(yvars) .* randn(length(xdata))

fit = curve_fit(model, xdata, ydata, 1 ./ yvars, [0.5, 0.5])
println("norm(fit.param - [1.0, 2.0]) < 0.05 ? ", norm(fit.param - [1.0, 2.0]))
@assert norm(fit.param - [1.0, 2.0]) < 0.05
@test norm(fit.param - [1.0, 2.0]) < 0.05
@test fit.converged

# can also get error estimates on the fit parameters
errors = margin_error(fit, 0.1)
println("norm(errors - [0.017, 0.075]) < 0.1 ?", norm(errors - [0.017, 0.075]))
@assert norm(errors - [0.017, 0.075]) < 0.1
@test norm(errors - [0.017, 0.075]) < 0.1

# test with user-supplied jacobian and weights
fit = curve_fit(model, jacobian_model, xdata, ydata, 1 ./ yvars, p0)
println("norm(fit.param - [1.0, 2.0]) < 0.05 ? ", norm(fit.param - [1.0, 2.0]))
@assert norm(fit.param - [1.0, 2.0]) < 0.05
@test norm(fit.param - [1.0, 2.0]) < 0.05
@test fit.converged

# Parameters can also be inferred using arbitrary precision
fit = curve_fit(model, xdata, ydata, 1 ./ yvars, BigFloat.(p0); x_tol=1e-20, g_tol=1e-20)
@test fit.converged
Expand Down
Loading