Global fit to multiple datasets

[Stéphane Gorsse]

Hello,

I have the code below that works well to optimize parameters (sigma_0, k1 and k2) to a set of given data by fitting a function that contains a numerically calculated integer.

 

# Constants

M = 3

b = 0.254e-9

alpha = 0.4

G = 83000

 

def f(rho, x, D,k1,k2):

    return M*(1/(b*D) + k1/b*(rho)**(1/2) - k2*rho)

 

def g(x, rho_0, D,k1,k2):

    """

    Solution to the ODE rho'(x) = f(rho,x,D,k1,k2) with initial condition rho(x=0) = rho_0

    """

    rho = odeint(f, rho_0, x, args=(D,k1,k2))

    return rho

 

def iso(x, D,k1,k2):

    return M*alpha*G*b*(g(x, rho_0, D,k1,k2))**(1/2)

 

def flow(x, sigma_0, D,k1,k2):

    return sigma_0 + iso(x,D,k1,k2)

 

# Create the dataset

D=10e-6

sigma_0 = 300

k1 = 0.015

k2 = 1.4

x = np.linspace(0, 1, 101)

data = flow(x, sigma_0, D,k1,k2)

data += 10*np.random.normal(size=data.shape)

 

# Define the parameters to be optimized

pars = mod.make_params()

pars['D'].set(10e-6, vary=False)

pars['sigma_0'].set(value=300, min=200, max=500)

pars['k1'].set(value=0.01, min=0.001, max=0.1)

pars['k2'].set(value=4, min=1, max=9)

 

# Run the global fit to all the data sets

mod = Model(flow)

result = mod.fit(data,pars,x=x)

print(result.fit_report())

 

fig = plt.figure(figsize=(6,5))

plt.scatter(x, data, color='gray', s=80,zorder=1)

plt.plot(x, result.best_fit.transpose().flatten(), 'r-', label='best fit')

plt.show()

 

Now I would like to adapt my code to perform a global fit to find the best common parameters (sigma_0, k1 and k2) that fit best to multiple datasets related to different values of parameter D. Question: could you suggest an approach to achieve this goal?

 

# Create the dataset1

D_1=10e-6

sigma_0 = 300

k1 = 0.015

k2 = 1.4

x = np.linspace(0, 1, 101)

data_1 = flow(x, sigma_0, D_1,k1,k2)

data_1 += 10*np.random.normal(size=data.shape)

 

# Create the dataset2

D_2=50e-6

x = np.linspace(0, 1, 101)

data_2 = flow(x, sigma_0, D_2,k1,k2)

data_2 += 10*np.random.normal(size=data.shape)

 

# Run the global fit to all the data sets

# model1 = flow(x_global, params['sigma_0'], params['D_1'], params['k1'], params['k2'])

# model2 = flow(x_global, params['sigma_0'], params['D_2'], params['k1'], params['k2'])

 

fig = plt.figure(figsize=(6,5))

plt.scatter(x, data_1, color='gray', s=80)

plt.scatter(x, data_2, color='black', s=80)

#plt.plot(x, result.best_fit.transpose().flatten(), 'r-', label='best fit')

plt.show()

Thanks in advance for your help

Stéphane

[Renee Otten]

Hi Stephane, 

did you read the documentation? To me it seems you want to follow the approach described here:

https://lmfit.github.io/lmfit-py/faq.html#how-can-i-fit-multiple-data-sets 

and an example is given here:

https://lmfit.github.io/lmfit-py/examples/example_fit_multi_datasets.html

Take a look at that first please and try to adapt it to your model; if you run into specific questions, please post the code you have here so that we can take a look.

Good luck,

Renee

--
You received this message because you are subscribed to the Google Groups "lmfit-py" group.
To unsubscribe from this group and stop receiving emails from it, send an email to lmfit-py+u...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/lmfit-py/f5b3cf63-b922-4561-8261-c28859e5e2e9n%40googlegroups.com.

[Renee Otten]

hi again Stephane,

I was looking at your question again and see that you are using the Model class. Just to add to my earlier response; if you look at the documentation/example you’ll see that these use the Minimizer class directly. 

I don’t think you can easily do the type of global fit using the Model class; but Matt can correct me if I’m wrong about this. So I’d suggest re-writing your code to use the Minimizer class and then follow the approach shown in the links I mentioned before.  

Give it a try and if you have trouble converting to use the Minimizer class, let us know. 

Renee

On Aug 7, 2021, at 2:26 PM, Renee Otten <otten...@gmail.com> wrote:

Hi Stephane, 

did you read the documentation? To me it seems you want to follow the approach described here:

https://lmfit.github.io/lmfit-py/faq.html#how-can-i-fit-multiple-data-sets 

and an example is given here:

[Mark Dean]

Hi,

I think you can use the model interface as follows, although one does have to construct a function.

import matplotlib.pyplot as plt
import numpy as np

import lmfit
from lmfit.lineshapes import gaussian

def multigauss(x, a0=1, a1=1, a2=1, a3=2, a4=3,
               c0=.1, c1=.2, c2=.4, c3=.6, c4=.1,
               s=.2):
    M = np.array([gaussian(x, a0, c0, s),
                  gaussian(x, a1, c1, s),
                  gaussian(x, a2, c2, s),
                  gaussian(x, a3, c3, s),
                  gaussian(x, a4, c4, s)])
    return M.T


x = np.linspace(-1, 2, 51)
data = []
for _ in np.arange(5):
    amp = 0.60 + 9.50*np.random.rand()
    cen = -0.20 + 1.20*np.random.rand()
    sig = 0.25 + 0.03*np.random.rand()
    dat = gaussian(x, amp, cen, sig) + np.random.normal(size=x.size, scale=0.1)
    data.append(dat)
data = np.array(data).T

model = lmfit.Model(multigauss)

params = model.make_params()

for k in params.keys():
    params[k].set(min=-5, max=10)


result = model.fit(data, x=x, params=params)
result

yfits = multigauss(x, **result.best_values)

plt.figure()

colors = plt.matplotlib.cm.nipy_spectral(np.linspace(0, 1, data.shape[1]))
for data_set, fit_set, color in zip(data.T, yfits.T, colors):
    plt.plot(x, fit_set, '-', color=color)
    plt.plot(x, data_set, '.', color=color)
 
plt.show()

To view this discussion on the web visit https://groups.google.com/d/msgid/lmfit-py/B82C671A-9748-49C0-86E6-085ECCF56BD5%40gmail.com.

--

--------------------------------------------

Mark P. M. Dean

Brookhaven National Laboratory

Condensed Matter Physics Bldg 734

Upton NY 11973

USA

Phone: 001 631 344 7847

Email: md...@bnl.gov

[Stéphane Gorsse]

Dear Renee, Mark,

Thank you very much for your advice. I modified my code to use the Minimizer class. However, I experience one trouble with the code below:

TypeError: type object got multiple values for keyword argument 'fcn_args'

I tried various things but I am stuck.

# Constants

M = 3

b = 0.254e-9

alpha = 0.4

G = 83000

def f(rho, x, D,k1,k2):

    return M*(1/(b*D) + k1/b*(rho)**(1/2) - k2*rho)

def g(x, rho_0, D,k1,k2):

    """

    Solution to the ODE rho'(x) = f(rho,x,k1,k2) with initial condition rho(x=0) = rho_0

    """

    rho = odeint(f, rho_0, x, args=(D,k1,k2))

    return rho

def iso(x, D,k1,k2):

    return M*alpha*G*b*(g(x, rho_0, D,k1,k2))**(1/2)

def flow(x, sigma_0, D,k1,k2):

    return sigma_0 + iso(x,D,k1,k2)

# Create the dataset1

D_1=10e-6

sigma_0 = 300

k1 = 0.015

k2 = 1.4

x = np.linspace(0, 1, 101)

data_1 = flow(x, sigma_0, D_1,k1,k2)

data_1 += 10*np.random.normal(size=data.shape)

# Create the dataset2

D_2=50e-6

x = np.linspace(0, 1, 101)

data_2 = flow(x, sigma_0, D_2,k1,k2)

data_2 += 10*np.random.normal(size=data.shape)

datasets = np.concatenate((data_1, data_2))

params =  Parameters()

params.add('sigma_0', 260, vary=True)

params.add('D_1', 10e-6, vary=False)

params.add('D_2', 50e-6, vary=False)

params.add('k1', 0.015, vary=True)

params.add('k2', 1.4, vary=True)

#def residual(params, x=None, y_1150_global=None, y_1100_global=None):

def residual(params, x, datasets):

    model1 = flow(x, params['sigma_0'], params['D_1'], params['k1'], params['k2'])

    model2 = flow(x, params['sigma_0'], params['D_2'], params['k1'], params['k2'])

    resid1 = data_1 - model1

    resid2 = data_2 - model2

    return np.concatenate((resid1, resid2))

fit = lmfit.minimize(residual, params, fcn_args=(x, datasets))

print(result.fit_report())                           

fig = plt.figure(figsize=(6,5))

plt.scatter(x, data_1, color='gray', s=80)

plt.scatter(x, data_2, color='black', s=80)

plt.show()

Cheers

Stéphane

[Renee Otten]

Hi Stephane, 

Unfortunately the code you posted doesn’t seem to be a “minimal, working example” in the sense that I cannot run it here. The “import numpy as np” is easily fixed, but then I get "NameError: name 'rho_0' is not defined”.

Just looking at your code, it seems to me that you probably don’t want to use “concatenate” on your data and then also the definition of your “residual” function should change. Perhaps the code below would work, but as I said I cannot test it since your example doesn’t run….

datasets = [data_1, data_2]

def residual(params, x, datasets):

    model1 = flow(x, params['sigma_0'], params['D_1'], params['k1'], params['k2'])

    model2 = flow(x, params['sigma_0'], params['D_2'], params['k1'], params['k2'])

    resid1 = datasets[0] - model1

    resid2 = datasets[1] - model2

    return np.concatenate((resid1, resid2))

Best, 

Renee

To view this discussion on the web visit https://groups.google.com/d/msgid/lmfit-py/bf2fd455-ea97-468f-bf3e-f2857b87e279n%40googlegroups.com.

[Stéphane Gorsse]

Hi Renee, 

Sorry for the mistake and thank you for your patience.

Below is the corrected code that you can run:

import numpy as np

import matplotlib.pyplot as plt

from scipy.integrate import odeint

import lmfit as lmfit

from lmfit import minimize, Parameters, Parameter, report_fit, Model

# Constants

M = 3

b = 0.254e-9

alpha = 0.4

G = 83000

rho_0 = 1e11

x = np.linspace(0, 1, 101)

def f(rho, x, D,k1,k2):

    return M*(1/(b*D) + k1/b*(rho)**(1/2) - k2*rho)

def g(x, rho_0, D,k1,k2):

    """

    Solution to the ODE rho'(x) = f(rho,x,k1,k2) with initial condition rho(x=0) = rho_0

    """

    rho = odeint(f, rho_0, x, args=(D,k1,k2))

    return rho

def iso(x, D,k1,k2):

    return M*alpha*G*b*(g(x, rho_0, D,k1,k2))**(1/2)

def flow(x, sigma_0, D,k1,k2):

    return sigma_0 + iso(x,D,k1,k2)

# Create the dataset1

D_1=10e-6

sigma_0 = 300

k1 = 0.015

k2 = 1.4

data_1 = flow(x, sigma_0, D_1,k1,k2)

data_1 += 10*np.random.normal(size=data_1.shape)

# Create the dataset2

D_2=50e-6

data_2 = flow(x, sigma_0, D_2,k1,k2)

data_2 += 10*np.random.normal(size=data_2.shape)

datasets = [data_1, data_2]

params =  Parameters()

params.add('sigma_0', 260, vary=True)

params.add('D_1', 10e-6, vary=False)

params.add('D_2', 50e-6, vary=False)

params.add('k1', 0.015, vary=True)

params.add('k2', 1.4, vary=True)

#def residual(params, x=None, y_1150_global=None, y_1100_global=None):

def residual(params, x, datasets):

    model1 = flow(x, params['sigma_0'], params['D_1'], params['k1'], params['k2'])

    model2 = flow(x, params['sigma_0'], params['D_2'], params['k1'], params['k2'])

    resid1 = data_1 - model1

    resid2 = data_2 - model2

    return np.concatenate((resid1, resid2))

fit = lmfit.minimize(residual, params, fcn_args=(x, datasets))

print(result.fit_report())                           

fig = plt.figure(figsize=(6,5))

plt.scatter(x, data_1, color='gray', s=80)

plt.scatter(x, data_2, color='black', s=80)

plt.show()

Cheers

Stéphane

Le lundi 9 août 2021 à 13:37:45 UTC+2, Renee Otten a écrit :

Hi Stephane, 

Unfortunately the code you posted doesn’t seem to be a “minimal, working example” in the sense that I cannot run it here. The “import numpy as np” is easily fixed, but then I get "NameError: name 'rho_0' is not defined”.

Just looking at your code, it seems to me that you probably don’t want to use “concatenate” on your data and then also the definition of your “residual” function should change. Perhaps the code below would work, but as I said I cannot test it since your example doesn’t run….

datasets = [data_1, data_2]

def residual(params, x, datasets):

    model1 = flow(x, params['sigma_0'], params['D_1'], params['k1'], params['k2'])

    model2 = flow(x, params['sigma_0'], params['D_2'], params['k1'], params['k2'])

    resid1 = datasets[0] - model1

    resid2 = datasets[1] - model2

    return np.concatenate((resid1, resid2))

Best, 

Renee

[Renee Otten]

Hi Stephane, 

Okay, that was actually almost correct ;) 

The “minimize” function wants to have the arguments to the residual function as “args” (see: https://lmfit.github.io/lmfit-py/fitting.html#the-minimize-function). You were supplying them as “fcn_args”; the latter is the correct name when using the Minimizer class (see: https://lmfit.github.io/lmfit-py/fitting.html#using-the-minimizer-class).

I changed a few other minor things (it can likely  be cleaned up a bit more, but I will leave that for you); anywhon, the code below works for me and appears to give a good fit:

Hope that helps!

Renee

*****

import numpy as np

import matplotlib.pyplot as plt

from scipy.integrate import odeint

from lmfit import minimize, Parameters, report_fit

# Constants

M = 3

b = 0.254e-9

alpha = 0.4

G = 83000

rho_0 = 1e11

x = np.linspace(0, 1, 101)

def f(rho, x, D, k1, k2):

    return M * (1 / (b * D) + k1 / b * (rho) ** (1 / 2) - k2 * rho)

def g(x, rho_0, D, k1, k2):

    """

    Solution to the ODE rho'(x) = f(rho,x,k1,k2) with initial condition rho(x=0) = rho_0

    """

    rho = odeint(f, rho_0, x, args=(D, k1, k2))

    return rho

def iso(x, D, k1, k2):

    return M * alpha * G * b * (g(x, rho_0, D, k1, k2)) ** (1 / 2)

def flow(x, sigma_0, D, k1, k2):

    return np.ravel(sigma_0 + iso(x, D, k1, k2))

# Create the dataset1

D_1 = 10e-6

sigma_0 = 300

k1 = 0.015

k2 = 1.4

data_1 = flow(x, sigma_0, D_1, k1, k2)

data_1 += 10 * np.random.normal(size=data_1.shape)

# Create the dataset2

D_2 = 50e-6

data_2 = flow(x, sigma_0, D_2, k1, k2)

data_2 += 10 * np.random.normal(size=data_2.shape)

datasets = [data_1, data_2]

params = Parameters()

params.add("sigma_0", 260, vary=True)

params.add("D_1", 10e-6, vary=False)

params.add("D_2", 50e-6, vary=False)

params.add("k1", 0.015, vary=True)

params.add("k2", 1.4, vary=True)

def residual(params, x, datasets=None):

    model1 = flow(x, params["sigma_0"], params["D_1"], params["k1"], params["k2"])

    model2 = flow(x, params["sigma_0"], params["D_2"], params["k1"], params["k2"])

    if datasets is None:

        return (model1, model2)

    resid1 = datasets[0] - model1

    resid2 = datasets[1] - model2

    return np.concatenate((resid1, resid2))

result = minimize(residual, params, args=(x, datasets))

report_fit(result)

fig = plt.figure(figsize=(6, 5))

plt.scatter(x, data_1, color="gray", s=80)

plt.scatter(x, data_2, color="black", s=80)

fitted_lines = residual(result.params, x)

plt.plot(x, fitted_lines[0], label="fit dataset 1")

plt.plot(x, fitted_lines[1], label="fit dataset 2")

plt.legend()

plt.show()

[Stéphane Gorsse]

Hi Renee,

This is great. Thank you so much, your precious help is very much appreciated.

Cheers

Stéphane