Getting started with CurveFit.jl

This tutorial introduces the basic workflow of CurveFit by walking through several simple examples.

Linear Fitting

Fit a linear function y = a * x + b:

using CurveFit

# Generate sample data: y = 2.5 * x + 3.0
x = collect(0:0.1:10)
y = @. 2.5 * x + 3.0

# Create the problem and solve
prob = CurveFitProblem(x, y)
sol = solve(prob, LinearCurveFitAlgorithm())

# Access the coefficients: sol.u = (a, b)
println("Slope (a): ", sol.u[1])
println("Intercept (b): ", sol.u[2])

# Evaluate the solution at a point
println("Prediction at x=5: ", sol(5.0))

Slope (a): 2.500000000000001
Intercept (b): 2.999999999999998
Prediction at x=5: 15.500000000000002

Nonlinear Curve Fitting

For arbitrary nonlinear functions, use NonlinearCurveFitProblem:

using CurveFit

# Define a nonlinear function: y = a[1] + a[2] * x^a[3]
fn(a, x) = @. a[1] + a[2] * x^a[3]

# True parameters
true_params = [3.0, 2.0, 0.7]

# Generate sample data
x = collect(1.0:0.5:10.0)
y = fn(true_params, x)

# Create problem with initial guess for parameters
u0 = [0.5, 0.5, 0.5]
prob = NonlinearCurveFitProblem(fn, u0, x, y)
sol = solve(prob)

println("Fitted parameters: ", sol.u)
println("Prediction at x=5: ", sol(5.0))

Fitted parameters: [3.0, 2.0, 0.7]
Prediction at x=5: 9.170338627200096

Using ScalarModel for Scalar Functions

The @. macro in the example above broadcasts the function over all data points. If you prefer to define a simpler scalar function that operates on one data point at a time, use ScalarModel:

using CurveFit

# Define a scalar function (no @. needed): y = a[1] + a[2] * x^a[3]
fn_scalar(a, x) = a[1] + a[2] * x^a[3]

# True parameters
true_params = [3.0, 2.0, 0.7]

# Generate sample data
x = collect(1.0:0.5:10.0)
y = fn_scalar.(Ref(true_params), x)  # Note: we still need to broadcast for data generation

# Create problem with ScalarModel wrapper
u0 = [0.5, 0.5, 0.5]
prob = NonlinearCurveFitProblem(ScalarModel(fn_scalar), u0, x, y)
sol = solve(prob)

println("Fitted parameters: ", sol.u)
println("Prediction at x=5: ", sol(5.0))

Fitted parameters: [3.0, 2.0, 0.7]
Prediction at x=5: 9.170338627200096

Migration from LsqFit.jl

If you're migrating from LsqFit.jl, note these key differences:

Parameter order is reversed: LsqFit uses model(x, p), CurveFit uses model(p, x).
Vectorized functions: CurveFit expects functions that operate on array data. Using ScalarModel makes migration easier.

Polynomial Fitting

Fit a polynomial of a given degree:

using CurveFit

# Generate sample data: y = 1.0 + 2.0*x + 3.0*x^2
x = collect(range(1, stop=10, length=20))
y = @. 1.0 + 2.0 * x + 3.0 * x^2

# Create the problem and solve with degree 2 polynomial
prob = CurveFitProblem(x, y)
sol = solve(prob, PolynomialFitAlgorithm(degree=2))

# Access the coefficients: [c0, c1, c2] for c0 + c1*x + c2*x^2
println("Coefficients: ", sol.u)

# Evaluate the solution at a point
println("Prediction at x=5: ", sol(5.0))

Coefficients: [1.0000000000000715, 1.999999999999975, 3.0000000000000018]
Prediction at x=5: 85.99999999999999

Exponential Fitting

Fit an exponential function y = b * exp(a * x):

using CurveFit

# Generate sample data: y = 2.0 * exp(0.3 * x)
x = collect(range(0, stop=5, length=20))
y = @. 2.0 * exp(0.3 * x)

prob = CurveFitProblem(x, y)
sol = solve(prob, ExpCurveFitAlgorithm())

# sol.u[1] = a (exponent coefficient)
# sol.u[2] = log(b)
println("Exponent coefficient (a): ", sol.u[1])
println("Scale factor (b): ", exp(sol.u[2]))

Exponent coefficient (a): 0.2999999999999999
Scale factor (b): 2.000000000000001

Power Law Fitting

Fit a power law y = b * x^a:

using CurveFit

# Generate sample data: y = 2.0 * x^0.8
x = collect(range(1, stop=10, length=20))
y = @. 2.0 * x^0.8

prob = CurveFitProblem(x, y)
sol = solve(prob, PowerCurveFitAlgorithm())

# sol.u[1] = a (exponent)
# sol.u[2] = log(b)
println("Exponent (a): ", sol.u[1])
println("Scale factor (b): ", exp(sol.u[2]))

Exponent (a): 0.8000000000000014
Scale factor (b): 1.999999999999995

Sum of Exponentials

Fit a sum of exponentials: y = k + p * exp(λ * t):

using CurveFit

# Generate sample data: y = 2.0 + 3.0*exp(-0.5*t)
t = collect(range(0, stop=10, length=50))
y = @. 2.0 + 3.0 * exp(-0.5 * t)

prob = CurveFitProblem(t, y)
sol = solve(prob, ExpSumFitAlgorithm(n=1, withconst=true))

# Access fitted parameters (returned as arrays for n exponentials)
# k[] extracts the scalar from a 1-element array
println("Constant (k): ", sol.u.k[])
println("Amplitude (p): ", sol.u.p[])
println("Decay rate (λ): ", sol.u.λ[])

Constant (k): 1.9995923847296964
Amplitude (p): 2.999598865560582
Decay rate (λ): -0.499566604252092

Modified King fitting

Fit with the modified king law: x^2 = a + b * y^n

using CurveFit

# Generate the data: x^2 = a + b * y^n
x = collect(10.0:20.0)
θ_ref = [5.0, 1.3, 2.5]
y = @. ((x^2 - θ_ref[2])/θ_ref[3])^(1/θ_ref[1])

# Works with and without an initial guess
prob = CurveFitProblem(x, y)
sol = solve(prob, ModifiedKingCurveFitAlgorithm())

println("Fitted parameters: ", sol.u)

Fitted parameters: [1.2999999999996883, 2.5000000000000173, 4.999999999999994]

Rational polynomial fitting

Fit a rational function: y = p(x)/q(x)

using CurveFit

# Generate sample data: y = (1 + 2*x) / (1 + 0.5*x - 0.1*x^2)
x = collect(range(0, stop=5, length=30))
y = @. (1.0 + 2.0*x) / (1.0 + 0.5*x - 0.1*x^2)

# Create the problem and solve with numerator degree 1, denominator degree 2
prob = CurveFitProblem(x, y)
alg = RationalPolynomialFitAlgorithm(num_degree=1, den_degree=2)
sol = solve(prob, alg)

# Access fitted coefficients
# Numerator: p0 + p1*x
println("Numerator coefficients: ", sol.u[1:2])
# Denominator: 1 + q1*x + q2*x^2
println("Denominator coefficients: ", vcat(1.0, sol.u[3:4]))

# Evaluate the solution at a point
println("Prediction at x=2: ", sol(2.0))

Numerator coefficients: [1.0, 2.0000000000000004]
Denominator coefficients: [1.0, 0.5000000000000001, -0.10000000000000002]
Prediction at x=2: 3.1250000000000004