2
$\begingroup$

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.

Improve this question
$\endgroup$
13
  • 1
    $\begingroup$ You don't optimize models, you optimize an objective (cost) function. What is your objective function? $\endgroup$
    – Roland
    May 29, 2018 at 14:30
  • $\begingroup$ 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. $\endgroup$
    – James Phillips
    May 29, 2018 at 19:11
  • $\begingroup$ @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. $\endgroup$
    – user9862881
    May 29, 2018 at 19:27
  • 1
    $\begingroup$ 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. $\endgroup$
    – Roland
    May 30, 2018 at 5:40
  • 1
    $\begingroup$ @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. $\endgroup$
    – user9862881
    May 30, 2018 at 11:12

1 Answer 1

Reset to default
2
$\begingroup$

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)
Improve this answer
$\endgroup$
2
  • 1
    $\begingroup$ thank you for the code. I will see if I can adapt it into r. Alternately I should probably learn Python at some point. $\endgroup$
    – user9862881
    May 30, 2018 at 7:00
  • $\begingroup$ I found this response helpful and wanted to share that the key principle shown by this answer (appending two obj function outputs together as a single obj function) works for python's lmfit package as well. $\endgroup$
    – vector07
    Oct 22, 2020 at 23:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.