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.

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

Allow curve_fit() to take a data argument to be passed to model function #14

Closed
wants to merge 2 commits into from
Closed
Show file tree
Hide file tree
Changes from all commits
Commits
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
7 changes: 4 additions & 3 deletions README.md
Original file line number Diff line number Diff line change
Expand Up @@ -41,13 +41,14 @@ There are top-level methods `curve_fit()` and `estimate_errors()` that are usefu
Existing Functionality
----------------------

`fit = curve_fit(model, x, y, w, p0; kwargs...)`:
`fit = curve_fit(model, x, y, w, p0; margs=[], kwargs...)`:

* `model`: function that takes two arguments (x, params)
* `x`: the independent variable
* `y`: the dependent variable that constrains `model`
* `w`: weight applied to the residual; can be a vector (of `length(x)` size) or matrix (inverse covariance)
* `p0`: initial guess of the model parameters
* `margs`: constant array of arguments to be passed to the model function when set, model(x, params, margs)
* `kwargs`: tuning parameters for fitting, passed to `levenberg_marquardt` of `Optim.jl`
* `fit`: composite type of results (`LsqFitResult`)

Expand All @@ -62,7 +63,7 @@ This performs a fit using a non-linear iteration to minimize the (weighted) resi
* `alpha`: confidence limit to calculate for the errors on parameters
* `sigma`: typical (symmetric) standard deviation for each parameter

This returns the error or uncertainty of each parameter fit to the model and already scaled by the associated degrees of freedom. Please note, this is a LOCAL quantity calculated from the jacobian of the model evaluated at the best fit point and NOT the result of a parameter exploration.
This returns the error or uncertainty of each parameter fit to the model and already scaled by the associated degrees of freedom. Please note, this is a LOCAL quantity calculated from the jacobian of the model evaluated at the best fit point and NOT the result of a parameter exploration.

----

Expand All @@ -71,4 +72,4 @@ This returns the error or uncertainty of each parameter fit to the model and alr
* `fit`: result of curve_fit (a `LsqFitResult` type)
* `covar`: parameter covariance matrix calculated from the jacobian of the model at the fit point

This returns the parameter covariance matrix evaluted at the best fit point.
This returns the parameter covariance matrix evaluted at the best fit point.
26 changes: 20 additions & 6 deletions src/curve_fit.jl
Original file line number Diff line number Diff line change
Expand Up @@ -33,20 +33,30 @@ function lmfit(f::Function, p0; kwargs...)
return LsqFitResult(dof, p, f(p), g(p))
end

function curve_fit(model::Function, xpts, ydata, p0; kwargs...)
function curve_fit{T<:Any}(model::Function, xpts, ydata, p0; margs::Array{T,1}=[], kwargs...)
# construct the cost function
f(p) = model(xpts, p) - ydata
if length(margs)==0
f(p) = model(xpts, p) - ydata
else
f(p) = model(xpts, p, margs) - ydata
end

lmfit(f,p0; kwargs...)
end

function curve_fit(model::Function, xpts, ydata, wt::Vector, p0; kwargs...)
function curve_fit{T<:Any}(model::Function, xpts, ydata, wt::Vector, p0; margs::Array{T,1}=[], kwargs...)
# construct a weighted cost function, with a vector weight for each ydata
# for example, this might be wt = 1/sigma where sigma is some error term
f(p) = wt .* ( model(xpts, p) - ydata )
if length(margs)==0
f(p) = wt .* ( model(xpts, p) - ydata )
else
f(p) = wt .* ( model(xpts, p, margs) - ydata )
end

lmfit(f,p0; kwargs...)
end

function curve_fit(model::Function, xpts, ydata, wt::Matrix, p0; kwargs...)
function curve_fit{T<:Any}(model::Function, xpts, ydata, wt::Matrix, p0; margs::Array{T,1}=[], kwargs...)
# as before, construct a weighted cost function with where this
# method uses a matrix weight.
# for example: an inverse_covariance matrix
Expand All @@ -56,7 +66,11 @@ function curve_fit(model::Function, xpts, ydata, wt::Matrix, p0; kwargs...)
# This requires the matrix to be positive definite
u = chol(wt)

f(p) = u * ( model(xpts, p) - ydata )
if length(margs)==0
f(p) = u * ( model(xpts, p) - ydata )
else
f(p) = u * ( model(xpts, p, margs) - ydata )
end
lmfit(f,p0; kwargs...)
end

Expand Down