Simultaneously curve fitting for 2 models with shared parameters in R

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.

• You don't optimize models, you optimize an objective (cost) function. What is your objective function? – Roland May 29 '18 at 14:30
• I have done this using a single combined data set and if-else logic within a single combined model, but did not use R. for example, "if data_set_1 then model_1 else model_2" as pseudocode. In your case both y1 and y2 have the same number of data points, but if this is not the case later you might want to weight the data as a correction. – James Phillips May 29 '18 at 19:11
• @Roland I'm afraid that I may be using the wrong terminology, this is all quite new to me. I am looking to fit these models to my data by optimizing the shared parameters a and h. I have edited the question to make it clearer. – user9862881 May 29 '18 at 19:27
• What do you want to minimize? Just the sum of both sums of squares? That would only make some sense if both y-variables have the same unit. You must define the objective function to solve this. – Roland May 30 '18 at 5:40
• @Roland I have to admit that your questions just highlight how much I still don't know. I'm very much self taught so there are some, not insignificant, gaps in my knowledge. I was going on the assumption that I would be minimizing both sum of squares, but the 2 y variables do have different units. I actually don't know how I would go about defining the objective function but I am interested in learning this. I would be interested if you could point me in the direction of some starter resources. – user9862881 May 30 '18 at 11:12

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)

• thank you for the code. I will see if I can adapt it into r. Alternately I should probably learn Python at some point. – user9862881 May 30 '18 at 7:00