pkg load interval ## Datasets (observations or known error bounds) global X X = vertcat (infsup (0.1, 0.2), ... infsup (0.2, 0.3), ... infsup (0.8, 0.9)); global Y Y = vertcat (infsup (-0.2, 1.5), ... infsup (-0.3, 1.9), ... infsup (1.2, 2.5)); ## This function is to be inverted function [fval, a1, a2, a3, b1, b2] = F (y, a1, a2, a3, b1, b2) global X global Y # Vectorized forward evaluation A = [a1; a2; a3]; lnX = log (X); A_lnX = A .* lnX; b2_X = b2 .* X; b1_b2_X = b1 + b2_X; fval = A_lnX .* b1_b2_X; # Backward evaluation and refinement of the parameters y = intersect (y, fval); # Contract the factors of the multiplication A_lnX = mulrev (b1_b2_X, y, A_lnX); b1_b2_X = mulrev (A_lnX, y, b1_b2_X); A = mulrev (lnX, A_lnX, A); a1 = A(1, :); a2 = A(2, :); a3 = A(3, :); # Contract the sum b2_X = intersect (b2_X, b1_b2_X - b1); b1 = intersect (b1, union (b1_b2_X - b2_X, [], 1)); # Contract b2 b2 = union (mulrev (X, b2_X, b2), [], 1); endfunction ## Initial search ranges A = vertcat (infsup (3.1, 3.2), ... infsup (3.3, 3.4), ... infsup (3.5, 3.6)); B = vertcat (infsup (-10, 10), ... infsup (-10, 10)); [parameters, paving] = fsolve(@F, [A; B], Y, struct ("Contract", true)); a1 = parameters(1) # [3.1, 3.2001] a2 = parameters(2) # [3.2999, 3.4] a3 = parameters(3) # [3.5, 3.6001] b1 = parameters(4) # [0.55553, 1.8562] b2 = parameters(5) # [-9.3644, -4.5616] plot (paving(4, :), paving(5, :)) xlabel ("b1") ylabel ("b2")