Joint-fit using curve_fit with multiple equations and inputs - Is it even possible?

Question

I asked similar questions in January and April that @Miłosz Wieczór and @Joe were kind enough to show interest in. Now, I am facing a similar but different problem because I need to do a joint-fit with multiple equations and inputs in order to get the universal solutions for two parameters fc and alpha. My code (which is based on the answers from the previous questions) is as follows:

import numpy as np
from numpy import linalg
import math
from scipy.optimize import curve_fit, least_squares, minimize

ya_exp  = np.array([1, 1.3, 1.7, 2.1, 2.7, 3.5, 4.5, 5.8, 7.5, 9.7, 12, 16, 21, 27, 34, 44, 57, 73, 94, 120, 156, 250000])

yb_exp  = np.array([1, 1.3, 1.7, 2.1, 2.7, 3.5, 4.5, 5.8, 7.5, 9.7, 12, 16, 21, 27, 34, 44, 57, 73, 94, 120, 156])

xa1_exp = np.array([4.68, 4.70, 4.71, 4.72, 4.74, 4.75, 4.76, 4.77, 4.79, 4.80,
                    4.82, 4.83, 4.85, 4.87, 4.89, 4.90, 4.96, 4.99, 5.02, 5.06,
                    5.11, 6.23])

xb1_exp = np.array([0.018, 0.023, 0.019, 0.023, 0.022, 0.023, 0.023, 0.023, 0.023, 0.024,
                    0.025, 0.028, 0.032, 0.033, 0.034, 0.037, 0.040, 0.043, 0.045, 0.047,
                    0.049])

xa2_exp = np.array([7.01, 7.03, 7.04, 7.04, 7.04, 7.10, 7.13, 7.13, 7.16, 7.14,
                    7.19, 7.18, 7.19, 7.22, 7.24, 7.28, 7.32, 7.35, 7.40, 7.45,
                    7.49, 10.1])

xb2_exp = np.array([0.008, 0.009, 0.008, 0.009, 0.008, 0.010, 0.010, 0.010, 0.011, 0.012,
                    0.016, 0.017, 0.020, 0.023, 0.027, 0.029, 0.036, 0.040, 0.043, 0.046,
                    0.052])

xa3_exp = np.array([5.67, 5.67, 5.68, 5.69, 5.72, 5.74, 5.74, 5.76, 5.76, 5.79,
                    5.81, 5.81, 5.83, 4.86, 5.89, 5.91, 5.96, 6.00, 6.04, 6.10,
                    6.14, 7.56])

xb3_exp = np.array([0.011, 0.011, 0.012, 0.011, 0.012, 0.012, 0.012, 0.012, 0.015, 0.017,
                    0.021, 0.026, 0.028, 0.030, 0.036, 0.039, 0.046, 0.050, 0.056, 0.059,
                    0.063])

xa1_zero = np.min(xa1_exp)
xa1_inf  = np.max(xa1_exp)

xa2_zero = np.min(xa2_exp)
xa2_inf  = np.min(xa2_exp)

xa3_zero = np.min(xa3_exp)
xa3_inf  = np.min(xa3_exp)

ig_fc    = 500
ig_alpha = 0.35

def CCXA(f_exp, fc, alpha):
    x  = np.log(f_exp/fc)
    R  = xa_zero + 1/2 * (xa_inf - xa_zero) * (1 + np.sinh((1 - alpha) * x) / (np.cosh((1 - alpha) * x) + np.sin(1/2 * alpha * math.pi)))
    I  = 1/2 * (xa_inf - xa_zero) * np.cos(alpha * math.pi / 2) / (np.cosh((1 - alpha) * x) + np.sin(alpha * math.pi / 2))
    RI = np.sqrt(R ** 2 + I ** 2)
    return RI

def CCXB(f_exp, fc, alpha):
    x  = np.log(f_exp/fc)
    R  = xa_zero + 1/2 * (xa_inf - xa_zero) * (1 + np.sinh((1 - alpha) * x) / (np.cosh((1 - alpha) * x) + np.sin(1/2 * alpha * math.pi)))
    I  = 1/2 * (xa_inf - xa_zero) * np.cos(alpha * math.pi / 2) / (np.cosh((1 - alpha) * x) + np.sin(alpha * math.pi / 2))
    iQ = I / R
    return iQ

poptXA, pcovXA = curve_fit(CCXA, ya_exp, xa_exp, p0=(ig_fc, ig_alpha))

poptXB, pcovXB = curve_fit(CCXB, yb_exp, xb_exp, p0=(ig_fc, ig_alpha))

def objective(e_exp, f_exp):
    poptXA, pcovXA = curve_fit(CCXA, ya_exp, xa_exp, p0=(ig_fc, ig_alpha))

    poptXB, pcovXB = curve_fit(CCXB, yb_exp, xb_exp, p0=(ig_fc, ig_alpha))
    
    err_total = np.sum(np.sqrt(np.diag(pcovXA))) + np.sum(np.sqrt(np.diag(pcovXB)))

    delta = linalg.norm(poptXB - poptXA)
    
    return err_total, delta

test = objective(xa_exp, ya_exp)

My first problem is that I am not sure how to make CCXA and CCXB search and locate xa_exp and xb_exp from the global scope because the different variables are defined by other names: xa1_exp, xa2_exp, and xa3_exp plus xb1_exp, xb2_exp, and xb2_exp. Again, I am interested in fc and alpha as optimizable parameters. curve_fit(CCXA, ya_exp, xa_exp, p0=(ig_fc, ig_alpha)) is able to optimize for fc and alpha but relies on xa_exp and xb_exp from the global scope. How can I pass these to the functions when they are different? Also, notice that the length of ya_exp, xa1_exp, xa2_exp, and xa3_exp is 22 while the length of yb_exp, xb1_exp, xb2_exp, and xb3_exp is 21.

My second problem is that I have no idea how to write a function that uses curve_fit as a global joint-fit that includes all the values. In other words, I want to find the best universal fits for fc and alpha for all the all the values, so that I get one global (or universal?) fit instead of six independent fits. objective provides test[0]=poptXB[0]-poptXA[0], while test[1] is the norm of popXA and poptXB, but neither returns provide the fitted values for fc and alpha, which is what I seek.

Is this possible?

EDIT 26.08.2020 (and 30.08.2020)

I came across another question regarding joint fit and adjusted my code accordingly:

from lmfit import minimize, Parameters, fit_report

params = Parameters()
params.add('fc', value=500)
params.add('alpha', value=0.2)

def CCXA(f_exp, xa_zero, xa_inf, fc, alpha):
    x  = np.log(f_exp/fc)
    R  = xa_zero + 1/2 * (xa_inf - xa_zero) * (1 + np.sinh((1 - alpha) * x) / (np.cosh((1 - alpha) * x) + np.sin(1/2 * alpha * math.pi)))
    I  = 1/2 * (xa_inf - xa_zero) * np.cos(alpha * math.pi / 2) / (np.cosh((1 - alpha) * x) + np.sin(alpha * math.pi / 2))
    RI = np.sqrt(R ** 2 + I ** 2)
    return RI

def CCXB(f_exp, xa_zero, xa_inf, fc, alpha):
    x  = np.log(f_exp/fc)
    R  = xa_zero + 1/2 * (xa_inf - xa_zero) * (1 + np.sinh((1 - alpha) * x) / (np.cosh((1 - alpha) * x) + np.sin(1/2 * alpha * math.pi)))
    I  = 1/2 * (xa_inf - xa_zero) * np.cos(alpha * math.pi / 2) / (np.cosh((1 - alpha) * x) + np.sin(alpha * math.pi / 2))
    iQ = I / R
    return iQ

def fit_function(params, ya_data=None, xa1_data=None, xa2_data=None, xa3_data=None, yb_data=None, xb1_data=None, xb2_data=None, xb3_data=None):
    xa_data = np.array([xa1_data, xa2_data, xa3_data])
    modxa   = np.array([CCXA(ya_data, np.min(i), np.max(i), params['fc'], params['alpha']) for i in xa_data])
    #sigma   = np.array([1 / np.sqrt(i+1) for i in xa_data])
    #sigma   = np.full((1,22), 0.5)
    #resxa   = np.array([(i - j)/k for i, j, k in zip(xa_data, modxa, sigma)])
    resxa   = np.array([i - j for i, j in zip(xa_data, modxa)])
    
    xb_data = np.array([xb1_data, xb2_data, xb3_data]) 
    modxb   = np.array([CCXB(yb_data, np.min(i), np.max(i), params['fc'], params['alpha']) for i in xa_data])
    #sigmb   = np.array([1 / np.sqrt(i+1) for i in xb_data])
    #sigmb   = np.full((1,21), 1)
    #resxb   = np.array([(i - j)/k for i, j, k in zip(xb_data, modxb, sigmb)])
    resxb   = np.array([i - j for i, j in zip(xb_data, modxb)])
    
    return np.concatenate((resxa.ravel(), resxb.ravel()))

Minor modifications of CCXA and CCXB plus the introduction of fit_function that is solved using lmfit provided values for fc and alpha:

[[Fit Statistics]]
    # fitting method   = leastsq
    # function evals   = 23
    # data points      = 129
    # variables        = 2
    chi-square         = 8.89790022
    reduced chi-square = 0.07006221
    Akaike info crit   = -340.945624
    Bayesian info crit = -335.225999
[[Variables]]
    fc:     1149.59953 +/- 572.031178 (49.76%) (init = 500)
    alpha:  0.64118507 +/- 0.04236375 (6.61%) (init = 0.5)
[[Correlations]] (unreported correlations are < 0.100)
    C(fc, alpha) =  0.874

Since I am new to lmfit, I would like to know if this is how it is supposed to be used. Is what I did correct or is it completely wrong?

score 1 · Answer 1 · answered Aug 25 '20 at 14:05

1

It might be easier to use least_squares directly. It takes a vector of residuals, where you simply specify the l.h.s. of your two equations

answered Aug 25 '20 at 14:05

ev-br

24,968
9
65
78

Thank you for your response, @ev-br. I am not sure that I follow. Since I have two equations, I also have two left-hand sides. How exactly do I combine them using `least_squares` from `scipy`? I apologize for my ignorance. – naughty_waves Aug 25 '20 at 19:05
The vector of residuals has the length of twice the number of observations – ev-br Aug 26 '20 at 07:37

score 1 · Accepted Answer · answered Aug 27 '20 at 12:16

1

yes, your use of lmfit and np.concatenate() looks right to me to find the single set of parameter values that will minimize the residual of (model1, dataset1) and (model2, dataset2). I'm not entirely sure that I understand your model and the 3 x values. I think it might be done more efficiently, but that's sort of a different question of whether the fitting mechanics are doing what you intend - I think they are.

The next (trickier) question is whether you want to emphasize a misfit in one dataset more than the other. As written, your function asserts that all data points are equally weighted - that is a fine default approach, but may not always be the case.

answered Aug 27 '20 at 12:16

M Newville

7,486
2
16
29

Thank you for clarifying, @M Newville. I am also sure it could be done more efficient but I am a novice. I thought about using weights as well to "prioritize" `xb_data` over `xa_data` in the joint-fit procedure. How would I do that with this many datasets? Do I have to put individual weights on each of the six datasets? – naughty_waves Aug 29 '20 at 10:28
yes, you could simply multiply the residuals for each dataset by a different weighting factor to make a misfit in one dataset more important than a misfit in another dataset. it is not uncommon to weight by `1/uncertainty_in_data` (if you have those values) so that a misfit in noisier or "less certain" data is less important than a misfit in cleaner, more certain data. – M Newville Aug 29 '20 at 12:25
Thank you again, @M Newville. Despite feeling bad for taking so much of your time, I have another question. Unlike `Model`, `minimize` does not take any weights directly as inputs. It also seems that my problem is better suited for `minimize` and not `Model`. I need to define some sort of `sigma` and return `resxa=(data-model)/sigma` instead of `resxa=data-model`. How do I implement this if I want to focus the fit on primarily the `xb_data` and secondary the `xa_data`? Do I, for example, define `sigma=np.full((1,21), 0.5)` and `sigmb=np.full((1,21), 1` for `resxa` and `resxb`? – naughty_waves Aug 30 '20 at 12:46
with `minimize` your objective function returns the array to be minimized, so you will have to apply any weighting yourself. So, concatenate the appropriate `(data-model)/sigma` values for the two different datasets. – M Newville Aug 30 '20 at 17:19

Joint-fit using curve_fit with multiple equations and inputs - Is it even possible?

2 Answers2