0

I am trying to fit different data set with different non-linear function that shared some parameters and it look something like this:

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

#some non-linear function
def Sigma1x(x,C11,C111,C1111,C11111):
    return C11*x+1/2*C111*pow(x,2)+1/6*C1111*pow(x,3)+1/24*C11111*pow(x,4)

def Sigma2x(x,C12,C112,C1112,C11112):
    return C12*x+1/2*C112*pow(x,2)+1/6*C1112*pow(x,3)+1/24*C11112*pow(x,4)

def Sigma1y(y,C12,C111,C222,C112,C1111,C1112,C2222,C12222):
    return C12*y+1/2*(C111-C222+C112)*pow(y,2)+1/12*(C111+2*C1112-C2222)*pow(y,3)+1/24*C12222*pow(y,4)

def Sigma2y(y,C11,C222,C222,C2222):
    return C11*y+1/2*C222*pow(y,2)+1/6*C2222*pow(y,3)+1/24*C22222*pow(y,4)

def Sigmaz(z,C11,C12,C111,C222,C112,C1111,C1112,C2222,C1122,C11111,C11112,C122222,C11122,C22222):
    return (C11+C12)*z+1/2*(2*C111-C222+3*C112)*pow(z,2)+1/6*(3/2*C1111+4*C1112-1/2*C222+3*C1122)*pow(z,3)+\
                    1/24*(3*C11111+10*C11112-5*C12222+10*C11122-2*C22222)*pow(z,4)

# Experimental datasets

Xdata=np.loadtxt('x-direction.txt') #This contain x axis and two other dataset, should be fitted with Sigma1x and Sigma2x
Ydata=np.loadtxt('y-direction.txt') #his contain yaxis and two other dataset, should be fitted with Sigma1yand Sigma2y
Zdata=nploadtxt('z-direction.txt')#This contain z axis and one dataset  fitted with Sigmaz

The question is how to use optimize.leastsq or other packages to fit the data with the appropriate function, knowing that they share multiple paramaters?

Zack_Az
  • 36
  • 6
  • maybe use a loop? anyway, your post is unreadable... – Rafaó May 25 '20 at 10:07
  • The idea is to be able to define a Global fitting function that incorporate all the fitting parameters... Something like this https://www.originlab.com/doc/Tutorials/Multi-Functions-Global-Fitting – Zack_Az May 25 '20 at 12:50
  • please provide any code you have tried. This site is for coding problems & not a general advice forum. – DrBwts May 25 '20 at 14:59
  • @DrBwts I have modified the question with the code, but as a new to fitting in python, I am just constructing pieces that fits my requirement. Thanks for your understanding :) – Zack_Az May 25 '20 at 22:48

1 Answers1

0

I was able to solve ( partially the initial question). I found symfit a very comprehensive and easy to use. So i wrote the following code

 import matplotlib.pyplot as plt

 from symfit import *
 import numpy as np
 from symfit.core.minimizers import DifferentialEvolution, BFGS
 Y_strain = np.genfromtxt('Y_strain.csv', delimiter=',')
 X_strain=np.genfromtxt('X_strain.csv', delimiter=',')
 xmax=max(X_strain[:,0])
 xmin=min(X_strain[:,0])
 xdata = np.linspace(xmin, xmax, 50)
 ymax=max(Y_strain[:,0])
 ymin=max(Y_strain[:,0])
 ydata=np.linspace(ymin, ymax, 50)

 x,y,Sigma1x,Sigma2x,Sigma1y,Sigma2y=  variables('x,y,Sigma1x,Sigma2x,Sigma1y,Sigma2y')
 C11,C111,C1111,C11111,C12,C112,C1112,C11112,C222,C2222,C12222,C22222 =  parameters('C11,C111,C1111,C11111,C12,C112,C1112,C11112,C222,C2222,C12222,C22222')

model =Model({
     Sigma1x:C11*x+1/2*C111*pow(x,2)+1/6*C1111*pow(x,3)+1/24*C11111*pow(x,4),
     Sigma2x:C12*x+1/2*C112*pow(x,2)+1/6*C1112*pow(x,3)+1/24*C11112*pow(x,4),
     #Sigma1y:C12*y+1/2*(C111-C222+C112)*pow(y,2)+1/12*(C111+2*C1112-C2222)*pow(y,3)+1/24*C12222*pow(y,4),
     #Sigma2y:C11*y+1/2*C222*pow(y,2)+1/6*C2222*pow(y,3)+1/24*C22222*pow(y,4),  
})
  fit = Fit(model, x=X_strain[:,0], Sigma1x=X_strain[:,1],Sigma2x=X_strain[:,2])
  fit_result = fit.execute()
  print(fit_result)
  plt.scatter(Y_strain[:,0],Y_strain[:,2])
  plt.scatter(Y_strain[:,0],Y_strain[:,1])
  plt.plot(xdata, model(x=xdata, **fit_result.params).Sigma1x)
  plt.plot(xdata, model(x=xdata, **fit_result.params).Sigma2x)

However, The resulting fit is very bad :

 Parameter Value        Standard Deviation
 C11       1.203919e+02 3.988977e+00
 C111      -6.541505e+02 5.643111e+01
 C1111     1.520749e+03 3.713742e+02
 C11111    -7.824107e+02 1.015887e+03
 C11112    4.451211e+03 1.015887e+03
 C1112     -1.435071e+03 3.713742e+02
 C112      9.207923e+01 5.643111e+01
 C12       3.272248e+01 3.988977e+00
 Status message         Desired error not necessarily achieved due to precision loss.
 Number of iterations   59
 Objective              <symfit.core.objectives.LeastSquares object at 0x000001CC00C0A508>
 Minimizer              <symfit.core.minimizers.BFGS object at 0x000001CC7F84A548>

 Goodness of fit qualifiers:
 chi_squared            6.230510793023184
 objective_value        3.115255396511592
 r_squared              0.991979767376565

Any idea's how to improve the fit?

Zack_Az
  • 36
  • 6
  • did you plot the result? Why do you think it is a bad fit? The r_squared value looks pretty good to me. You can also try giving different initial guesses for the parameters if you already have some intuition for the result to help the fit along. – tBuLi Jun 01 '20 at 10:24
  • When the fit is ploted vs data, you can see that the fitting curve is very badly following the measured data[This is the results of the fit](https://i.stack.imgur.com/hy0J2.png) – Zack_Az Jun 01 '20 at 15:20
  • Perhaps you can use a global minimizer instead if finding good initial guesses is not that easy. See this example from the docs. Read it all the way, the finale is worth it ;). https://symfit.readthedocs.io/en/stable/examples/ex_mexican_hat.html?highlight=DifferentialEvolution#Global-minimization:-Skewed-Mexican-hat – tBuLi Jun 02 '20 at 10:00
  • Thanks for the suggestions when I add a minimizer i get the following error : bounds should be a sequence containing real valued (min, max) pairs for each value in x I added this line tot he code : fit = Fit(model, x=X_strain[:,0], Sigma1x=X_strain[:,1],Sigma2x=X_strain[:,2],y=Y_strain[:,0],Sigma1y=Y_strain[:,1],Sigma2y=Y_strain[:,2],minimizer=[DifferentialEvolution, BFGS]) – Zack_Az Jun 03 '20 at 07:29
  • Differential evolution only works if you set bounds for the parameters, see the documentation. – tBuLi Jun 03 '20 at 12:42
  • The code is correct, It's just my dataset that can't be fitted simultaneously with sharing prameters. Thanls for your answers. Is there a way to make the fitting faster? any idea how to parallelize the sequences? – Zack_Az Jun 09 '20 at 20:37