GSL::MultiFit::FdfSolver.new(T, n, p)This creates an instance of the GSL::MultiFit::FdfSolver class of type T for n observations and p parameters. The type T is given by a Fixnum constant or a String,
GSL::MultiFit::LMSDER or "lmsder"GSL::MultiFit::LMDER or "lmder"For example, the following code creates an instance of a Levenberg-Marquardt solver for 100 data points and 3 parameters,
solver = MultiFit::FdfSolver.new(MultiFit::LMDER, 100, 3)
GSL::MultiFit::FdfSolver#set(f, x)GSL::MultiFit::Function_fdf class (see below). The initial guess of the parameters x is given by a GSL::Vector object.GSL::MultiFit::FdfSolver#nameGSL::MultiFit::FdfSolver#xGSL::MultiFit::FdfSolver#dxGSL::MultiFit::FdfSolver#fGSL::MultiFit::FdfSolver#JGSL::MultiFit::FdfSolver#jacobianGSL::MultiFit::FdfSolver#jacgsl_multifit_nlin.h)GSL::MultiFit::Function_fdf.new()GSL::MultiFit::Function_fdf.new(f, df, p)GSL::MultiFit::Function_fdf.new(f, df, fdf, p)Function_fdf class, to a function with p parameters, The first two or three arguments are Ruby Proc objects to evaluate the function to minimize and its derivative (Jacobian).GSL::MultiFit::Function_fdf#set_procs(f, df, p)GSL::MultiFit::Function_fdf#set_procs(f, df, fdf, p)GSL::MultiFit::Function_fdf#set_data(t, y)GSL::MultiFit::Function_fdf#set_data(t, y, sigma)GSL::MultiFit::FdfSolver#iterateposition.GSL::MultiFit::FdfSolver#positionGSL::Vector object.GSL::MultiFit.covar(J, epsrel)GSL::MultiFit.covar(J, epsrel, covar)This method uses the Jacobian matrix J to compute the covariance matrix of the best-fit parameters. If an existing matrix covar is given, it is overwritten, and if not, this method returns a new matrix. The parameter epsrel is used to remove linear-dependent columns when J is rank deficient.
The covariance matrix is given by,
covar = (J^T J)^{-1}
and is computed by QR decomposition of J with column-pivoting. Any columns of R which satisfy
|R_{kk}| <= epsrel |R_{11}|
are considered linearly-dependent and are excluded from the covariance matrix (the corresponding rows and columns of the covariance matrix are set to zero).
A minimization procedure should stop when one of the following conditions is true:
The handling of these conditions is under user control. The method below allows the user to test the current estimate of the best-fit parameters.
GSL::MultiFit::FdfSolver#test_delta(epsabs, epsrel)This method tests for the convergence of the sequence by comparing the last step with the absolute error epsabs and relative error (epsrel to the current position. The test returns GSL::SUCCESS if the following condition is achieved,
|dx_i| < epsabs + epsrel |x_i|
for each component of x and returns GSL::CONTINUE otherwise.
GSL::MultiFit::FdfSolver#test_gradient(g, epsabs)GSL::MultiFit::FdfSolver#test_gradient(epsabs)This function tests the residual gradient g against the absolute error bound epsabs. If g is not given, it is calculated internally. Mathematically, the gradient should be exactly zero at the minimum. The test returns GSL::SUCCESS if the following condition is achieved,
\sum_i |g_i| < epsabs
and returns GSL::CONTINUE otherwise. This criterion is suitable for situations where the precise location of the minimum, x, is unimportant provided a value can be found where the gradient is small enough.
GSL::MultiFit::FdfSolver#gradientGSL::MultiFit.test_delta(dx, x, epsabs, epsrel)GSL::MultiFit.test_gradient(g, epsabs)GSL::MultiFit.gradient(jac, f, g)GSL::MultiFit.covar(jac, epsrel)GSL::MultiFit.covar(jac, epsrel, covar)GSL::MultiFit module.GSL::MultiFit::FdfSolver.fit(x, y, type[, guess])GSL::MultiFit::FdfSolver.fit(x, w, y, type[, guess])This method uses FdfSolver with the LMSDER algorithm to fit the data [x, y] to a function of type type. The returned value is an array of 4 elements, [coef, err, chisq, dof], where coef is an array of the fitting coefficients, err contains errors in estimating coef, chisq is the chi-squared, and dof is the degree-of-freedom in the fitting which equals to (data length - number of fitting coefficients). The optional argument guess is an array of initial guess of the coefficients. The fitting type type is given by a String as follows.
"gaussian": Gaussian fit, y = y0 + A exp(-(x-x0)^2/2/var), coef = [y0, A, x0, var]"gaussian_2peaks": 2-peak Gaussian fit, y = y0 + A1 exp(-(x-x1)^2/2/var1) + A2 exp(-(x-x2)^2/2/var2), coef = [y0, A1, x1, var1, A2, x2, var2]"exp": Exponential fit, y = y0 + A exp(-b x), coef = [y0, A, b]"dblexp": Double exponential fit, y = y0 + A1 exp(-b1 x) + A2 exp(-b2 x), coef = [y0, A1, b1, A2, b2]"sin": Sinusoidal fit, y = y0 + A sin(f x + phi), coef = [y0, A, f, phi]"lor": Lorentzian peak fit, y = y0 + A/((x-x0)^2 + B), coef = [y0, A, x0, B]"hill": Hill's equation fit, y = y0 + (m - y0)/(1 + (xhalf/x)^r), coef = [y0, n, xhalf, r]"sigmoid": Sigmoid (Fermi-Dirac) function fit, y = y0 + m/(1 + exp((x0-x)/r)), coef = [y0, m, x0, r]"power": Power-law fit, y = y0 + A x^r, coef = [y0, A, r]"lognormal": Lognormal peak fit, y = y0 + A exp[ -(log(x/x0)/width)^2 ], coef = [y0, A, x0, width]See Linear fitting for linear and polynomical fittings.
The following example program fits a weighted exponential model with background to experimental data, Y = A exp(-lambda t) + b. The first part of the program sets up the functions procf and procdf to calculate the model and its Jacobian. The appropriate fitting function is given by,
f_i = ((A exp(-lambda t_i) + b) - y_i)/sigma_i
where we have chosen t_i = i. The Jacobian matrix jac is the derivative of these functions with respect to the three parameters (A, lambda, b). It is given by,
J_{ij} = d f_i / d x_j
where x_0 = A, x_1 = lambda and x_2 = b.
require("gsl")
include GSL::MultiFit
# x: Vector, list of the parameters to determine
# t, y, sigma: Vectors, observational data
# f: Vector, function to minimize
procf = Proc.new { |x, t, y, sigma, f|
a = x[0]
lambda = x[1]
b = x[2]
n = t.size
for i in 0...n do
yi = a*Math::exp(-lambda*t[i]) + b
f[i] = (yi - y[i])/sigma[i]
end
}
# jac: Matrix, Jacobian
procdf = Proc.new { |x, t, y, sigma, jac|
a = x[0]
lambda = x[1]
n = t.size
for i in 0...n do
ti = t[i]
si = sigma[i]
ei = Math::exp(-lambda*ti)
jac.set(i, 0, ei/si)
jac.set(i, 1, -ti*a*ei/si)
jac.set(i, 2, 1.0/si)
end
}
f = MultiFit::Function_fdf.new(procf, procdf, 2)
# Create data
r = Rng.new()
t = Vector.new(n)
y = Vector.new(n)
sigma = Vector.new(n)
for i in 0...n do
t[i] = i
y[i] = 1.0 + 5*Math::exp(-0.1*t[i]) + r.gaussian(0.1)
sigma[i] = 0.1
end
f.set_data(t, y, sigma)
x = Vector.new(1.0, 0.0, 0.0) # initial guess
solver = FdfSolver.new(FdfSolver::LMSDER, n, np)
solver.set(f, x)
iter = 0
solver.print_state(iter)
begin
iter += 1
status = solver.iterate
solver.print_state(iter)
status = solver.test_delta(1e-4, 1e-4)
end while status == GSL::CONTINUE and iter < 500
covar = solver.covar(0.0)
position = solver.position
chi2 = pow_2(solver.f.dnrm2)
dof = n - np
printf("A = %.5f +/- %.5f\n", position[0], Math::sqrt(chi2/dof*covar[0][0]))
printf("lambda = %.5f +/- %.5f\n", position[1], Math::sqrt(chi2/dof*covar[1][1]))
printf("b = %.5f +/- %.5f\n", position[2], Math::sqrt(chi2/dof*covar[2][2]))
#!/usr/bin/env ruby
require("gsl")
include MultiFit
N = 100
y0 = 1.0
A = 2.0
x0 = 3.0
w = 0.5
r = Rng.alloc
x = Vector.linspace(0.01, 10, N)
sig = 1
# Lognormal function with noise
y = y0 + A*Sf::exp(-pow_2(Sf::log(x/x0)/w)) + 0.1*Ran::gaussian(r, sig, N)
guess = [0, 3, 2, 1]
coef, err, chi2, dof = MultiFit::FdfSolver.fit(x, y, "lognormal", guess)
y0 = coef[0]
amp = coef[1]
x0 = coef[2]
w = coef[3]
graph(x, y, y0+amp*Sf::exp(-pow_2(Sf::log(x/x0)/w)))