Solved – Simultaneously curve fitting for 2 models with shared parameters in R

curve fittingnlsr

Hi I am trying to curve fit 2 models (Van Genuchten & Mualem) with shared parameters in r. The models are:

# Models
mod1 = y1 ~ r+(s-r)/(1+(a*-h)^n)^(1-1/n)  
mod2 = y2 ~ (Ks*(1-((a*-h)^((1-1/n)*n))*((1+(a*-h)*n)^n)^-(1-1/n))^2)/((1+ 
(a*-h)^n)*(1-1/n)*b)

This is a curve fitting so I have values for the following:

#Parameters: These are currently guesstimates but will ultimately be 
#derived from experimental data.
s = 0.4; r = 0.1; b = 0.5; Ks = 0.1

# values - These are re-sampled to reduce the number of points
y1 = c(0.4735295, 0.4729359, 0.4719321,  0.4702538,  0.4674984,  0.4631038,  
0.4563907, 0.4467292,  0.4338216,  0.4179355,  0.3998694,  0.380652,  
0.3612183,  0.3422378,  0.324116,  0.307056, 0.2911322,  0.2763443,  
0.2626525,  0.2499978)

y2 = c(0.048356046, 0.038888484, 0.031257167, 0.025107757, 0.02015411, 
0.016165172, 0.012954387, 0.010371147, 0.00829387, 0.006624425, 
0.005283616, 0.00420754, 0.00334464, 0.002653328, 0.002100068, 0.001657819, 
0.001304786, 0.001023404, 0.000799527, 0.000621763) 

h = c(-1.258925, -1.995262, -3.162278, -5.011872, -7.943282, -12.589254, 
-19.952623, -31.622777, -50.118723, -79.432823, -125.892541, -199.526231, 
-316.227766, -501.187234, -794.328235, -1258.925412, -1995.262315, 
-3162.27766, -5011.872336, -7943.282347)

I have used nlsLM to model the curves individually but I have been unable to find a method to model together with the shared parameters "a" and "n". I would appreciate any advice.

Best Answer

Per my comment, here is working Python code that fits two data sets to two straight lines with different slopes and a single shared offset parameter. This is not intended as a direct answer, but is here so that I can post formatted code.

import numpy as np
import matplotlib
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

y1 = np.array([ 16.00,  18.42,  20.84,  23.26,  25.68])
y2 = np.array([-20.00, -25.50, -31.00, -36.50, -42.00])
comboY = np.append(y1, y2)

h = np.array([5.0, 6.1, 7.2, 8.3, 9.4])
comboX = np.append(h, h)


def mod1(data, a, b, c): # not all parameters are used here
        return a * data + c


def mod2(data, a, b, c): # not all parameters are used here
        return b * data + c


def comboFunc(comboData, a, b, c):
    # single data set passed in, extract separate data
    extract1 = comboData[:len(y1)] # first data
    extract2 = comboData[len(y2):] # second data

    result1 = mod1(extract1, a, b, c)
    result2 = mod2(extract2, a, b, c)

    return np.append(result1, result2)


# some initial parameter values
initialParameters = np.array([1.0, 1.0, 1.0])

# curve fit the combined data to the combined function
fittedParameters, pcov = curve_fit(comboFunc, comboX, comboY, initialParameters)

# values for display of fitted function
a, b, c = fittedParameters

y_fit_1 = mod1(h, a, b, c) # first data set, first equation
y_fit_2 = mod2(h, a, b, c) # second data set, second equation

plt.plot(comboX, comboY, 'D') # plot the raw data
plt.plot(h, y_fit_1) # plot the equation using the fitted parameters
plt.plot(h, y_fit_2) # plot the equation using the fitted parameters
plt.show()

print(fittedParameters)
Related Question