Trusted, Bureaucrats
4,223
edits
Changes
m
run by{{collapsible script| type = Python script| title = The <code>sympy_test.py</code> script.| intro = Run with: ''' <code>python sympy_test.py'''</code>| lang = python| script =# Tutorial from: http://scipy-lectures.github.io/advanced/sympy.html from sympy import *
</source>}} [[Category:Analysis techniques]]
added Category:Analysis techniques using HotCat
__TOC__
== Calculate Jacobian and Hessian matrix in python sympy for exponential decay function ==
See also:
Consider for example installing [[Epd_canopy | Enthought Canopy]]
== Tutorial with function for weighted difference between function evaluation with fitted parameters and measured values. ==
Created by:
SBiNLab
# In contrast to other Computer Algebra Systems, in SymPy you have to declare symbolic variables explicitly:
x = Symbol('x')
y = Symbol('y')
# You can differentiate any SymPy expression using diff(func, var). Examples:print("First some tests")f = sin(x)print("Function is", f)print("df/dx is:", diff(f, x)) f = x**2print("Function is", f)print("df/dx is:", diff(f, x)) print("") # Define exponential decay.i0 = Symbol('i0')times = Symbol('times')r2eff = Symbol('r2eff') print("Function to calculate intensity decay at time[i].")exp_dec = i0 * exp( -times * r2eff)print("Function f(r2eff, i0) is:", simplify(exp_dec))print("df/d_i0 should be: exp(-r2eff*times)", simplify(diff(exp_dec, i0)) )print("df/d_r2eff should be: -i0*times*exp(-r2eff*times)", simplify(diff(exp_dec, r2eff)) ) print("") print("Function to calculate the weighted difference between function evaluation with fitted parameters and measured values.")# Function to calculate the func.# The standard deviation of the measured values.errors = Symbol('errors')# The measured values.values = Symbol('values')# back_calc, the function which return the values, according to parameters.back_calc = exp_decfunc = 1.0 / errors * (values - back_calc)print("The func function is:", func) print("") print("Now calculate the Jacobian. The partial derivative matrix.\n")print("Jacobian is m rows with function derivatives and n columns of parameters.") # If you want to test for symbolic equality, one way is to subtract one expression from the other and run it through functions like expand(), simplify(), and trigsimp() and see if the equation reduces to 0.# https://github.com/sympy/sympy/wiki/Faq - Why does SymPy say that two equal expressions are unequal? print("Derivative of func function for exponential decay, with respect to r2eff")d_func_d_r2eff = diff(func, r2eff)print("d_func_d_r2eff=", d_func_d_r2eff) print("") print("Derivative of func function for exponential decay, with respect to i0")d_func_d_i0 = diff(func, i0)print("d_func_d_i0=", d_func_d_i0) print("\n") print("""Form the Jacobian matrix by:------------------------------------------------------------------------------from numpy import array, transpose d_func_d_r2eff = %sd_func_d_i0 = %sjacobian_matrix = transpose(array( [d_func_d_r2eff , d_func_d_i0] ) )------------------------------------------------------------------------------""" % (d_func_d_r2eff, d_func_d_i0) ) print("") print("Now calculate the Hessian. The second-order partial derivatives.\n")print("See for example: http://maxima-online.org/articles/hessian.html")print("If all second partial derivatives of f exist, then the Hessian matrix of f is the matrix:")print("Hf(x)i,j = d**2 f(x) / ( dx_i * dx_j )")print("We need to do:")print("diff(f(r2eff, i0), r2eff, r2eff)")print("diff(f(r2eff, i0), r2eff, i0)")print("diff(f(r2eff, i0), i0, r2eff)")print("diff(f(r2eff, i0), i0, i0)")d2_func_d_r2eff_d_r2eff = diff(func, r2eff, r2eff)d2_func_d_r2eff_d_i0 = diff(func, r2eff, i0)d2_func_d_i0_d_r2eff = diff(func, i0, r2eff)d2_func_d_i0_d_i0 = diff(func, i0, i0)print("d2_func_d_r2eff_d_r2eff=", d2_func_d_r2eff_d_r2eff)print("d2_func_d_r2eff_d_i0=", d2_func_d_r2eff_d_i0)print("d2_func_d_i0_d_r2eff=", d2_func_d_i0_d_r2eff)print("d2_func_d_i0_d_i0=", d2_func_d_i0_d_i0) print("\n") print("""Form the Hessian matrix by:------------------------------------------------------------------------------from numpy import array, transpose d2_func_d_r2eff_d_r2eff = %sd2_func_d_r2eff_d_i0 = %sd2_func_d_i0_d_r2eff = %sd2_func_d_i0_d_i0 = %shessian_matrix = transpose(array( [d2_func_d_r2eff_d_r2eff, d2_func_d_r2eff_d_i0, d2_func_d_i0_d_r2eff, d2_func_d_i0_d_i0] ) )------------------------------------------------------------------------------""" % (d2_func_d_r2eff_d_r2eff,d2_func_d_r2eff_d_i0, d2_func_d_i0_d_r2eff, d2_func_d_i0_d_i0) ) from sympy import Function, hessian, pprintf = Function('f')(r2eff, i0)}} == Tutorial with quadratic chi2 function ==Created by: Troels Emtekær Linnet <br>PhD student <source langbr>Copenhagen University <br>SBiNLab {{collapsible script| type ="Python"script| title = The <code>sympy_test.py</code> script.| intro = Run with: <code>python sympy_test.py</code>| lang = python| script =
# Tutorial from: http://scipy-lectures.github.io/advanced/sympy.html
------------------------------------------------------------------------------
""" % (d2_chi2_d_r2eff_d_r2eff,d2_chi2_d_r2eff_d_i0, d2_chi2_d_i0_d_r2eff, d2_chi2_d_i0_d_i0) )