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?

Answer 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.

Answer 2

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