Hi Troels, Could you derive the chi-squared Jacobian? Maybe the Jacobian I have been using is not correct - this is the one required for the Levenberg-Marquardt optimisation algorithm. Because the chi-squared is squared, its derivative will have a factor of 2 out the front, just like the gradient:
http://www.nmr-relax.com/manual/chi_squared_gradient.html It might be useful to add a Jacobian section to this part of the manual with the equations. Cheers, Edward On 28 August 2014 15:14, <tlin...@nmr-relax.com> wrote: > Author: tlinnet > Date: Thu Aug 28 15:14:16 2014 > New Revision: 25379 > > URL: http://svn.gna.org/viewcvs/relax?rev=25379&view=rev > Log: > Modified systemtest test Relax_disp.test_estimate_r2eff_err_methods() to show > the difference between using the direct function Jacobian, or the chi2 > function Jacobian. > > Added also the functionality to the estimate R2eff module, to switch between > using the different Jacobians. > > The results show, that R2eff can be estimated better. > > ---------------------- > The results are: > > R1rho at 799.8 MHz, for offset=118.078 ppm and dispersion point 431.0. > r2eff=8.646/8.646 r2eff_err=0.0348/0.0692 i0=202664.191/202664.191 > i0_err=699.6443/712.4201 > > R1rho at 799.8 MHz, for offset=118.078 ppm and dispersion point 651.2. > r2eff=10.377/10.377 r2eff_err=0.0403/0.0810 i0=206049.558/206049.558 > i0_err=776.4215/782.1833 > > R1rho at 799.8 MHz, for offset=118.078 ppm and dispersion point 800.5. > r2eff=10.506/10.506 r2eff_err=0.0440/0.0853 i0=202586.332/202586.332 > i0_err=763.9678/758.7052 > > R1rho at 799.8 MHz, for offset=118.078 ppm and dispersion point 984.0. > r2eff=10.903/10.903 r2eff_err=0.0476/0.0922 i0=203455.021/203455.021 > i0_err=837.8779/828.7280 > > R1rho at 799.8 MHz, for offset=118.078 ppm and dispersion point 1341.1. > r2eff=10.684/10.684 r2eff_err=0.0446/0.0853 i0=218670.412/218670.412 > i0_err=850.0210/830.9558 > > R1rho at 799.8 MHz, for offset=118.078 ppm and dispersion point 1648.5. > r2eff=10.501/10.501 r2eff_err=0.0371/0.0742 i0=206502.512/206502.512 > i0_err=794.0523/772.9843 > > R1rho at 799.8 MHz, for offset=124.247 ppm and dispersion point 1341.1. > r2eff=11.118/11.118 r2eff_err=0.0413/0.0827 i0=216447.241/216447.241 > i0_err=784.6562/788.0384 > > R1rho at 799.8 MHz, for offset=130.416 ppm and dispersion point 800.5. > r2eff=7.866/7.866 r2eff_err=0.0347/0.0695 i0=211869.715/211869.715 > i0_err=749.2776/763.6930 > > R1rho at 799.8 MHz, for offset=130.416 ppm and dispersion point 1341.1. > r2eff=9.259/9.259 r2eff_err=0.0331/0.0661 i0=217703.151/217703.151 > i0_err=682.2137/685.5838 > > R1rho at 799.8 MHz, for offset=130.416 ppm and dispersion point 1648.5. > r2eff=9.565/9.565 r2eff_err=0.0373/0.0745 i0=211988.939/211988.939 > i0_err=839.0313/827.0373 > > R1rho at 799.8 MHz, for offset=142.754 ppm and dispersion point 800.5. > r2eff=3.240/3.240 r2eff_err=0.0127/0.0253 i0=214417.382/214417.382 > i0_err=595.8865/613.4378 > > R1rho at 799.8 MHz, for offset=142.754 ppm and dispersion point 1341.1. > r2eff=5.084/5.084 r2eff_err=0.0177/0.0352 i0=226358.691/226358.691 > i0_err=660.5314/655.7670 > > R1rho at 799.8 MHz, for offset=179.768 ppm and dispersion point 1341.1. > r2eff=2.208/2.208 r2eff_err=0.0091/0.0178 i0=228620.553/228620.553 > i0_err=564.8353/560.0873 > > R1rho at 799.8 MHz, for offset=241.459 ppm and dispersion point 1341.1. > r2eff=1.711/1.711 r2eff_err=0.0077/0.0155 i0=224087.486/224087.486 > i0_err=539.4300/546.4217 > > Fitting with minfx to: 52V @N > ----------------------------- > > min_algor='Newton', c_code=True, constraints=False, chi2_jacobian?=False > ------------------------------------------------------------------------ > > R1rho at 799.8 MHz, for offset=118.078 ppm and dispersion point 431.0, with 4 > time points. r2eff=8.646 r2eff_err=0.0692, i0=202664.2, i0_err=712.4201, > chi2=3.758. > R1rho at 799.8 MHz, for offset=118.078 ppm and dispersion point 651.2, with 5 > time points. r2eff=10.377 r2eff_err=0.0810, i0=206049.6, i0_err=782.1833, > chi2=27.291. > R1rho at 799.8 MHz, for offset=118.078 ppm and dispersion point 800.5, with 5 > time points. r2eff=10.506 r2eff_err=0.0853, i0=202586.3, i0_err=758.7052, > chi2=13.357. > R1rho at 799.8 MHz, for offset=118.078 ppm and dispersion point 984.0, with 5 > time points. r2eff=10.903 r2eff_err=0.0922, i0=203455.0, i0_err=828.7280, > chi2=33.632. > R1rho at 799.8 MHz, for offset=118.078 ppm and dispersion point 1341.1, with > 5 time points. r2eff=10.684 r2eff_err=0.0853, i0=218670.4, i0_err=830.9558, > chi2=35.818. > R1rho at 799.8 MHz, for offset=118.078 ppm and dispersion point 1648.5, with > 5 time points. r2eff=10.501 r2eff_err=0.0742, i0=206502.5, i0_err=772.9843, > chi2=7.356. > R1rho at 799.8 MHz, for offset=124.247 ppm and dispersion point 1341.1, with > 5 time points. r2eff=11.118 r2eff_err=0.0827, i0=216447.2, i0_err=788.0384, > chi2=15.587. > R1rho at 799.8 MHz, for offset=130.416 ppm and dispersion point 800.5, with 5 > time points. r2eff=7.866 r2eff_err=0.0695, i0=211869.7, i0_err=763.6930, > chi2=14.585. > R1rho at 799.8 MHz, for offset=130.416 ppm and dispersion point 1341.1, with > 5 time points. r2eff=9.259 r2eff_err=0.0661, i0=217703.2, i0_err=685.5838, > chi2=79.498. > R1rho at 799.8 MHz, for offset=130.416 ppm and dispersion point 1648.5, with > 5 time points. r2eff=9.565 r2eff_err=0.0745, i0=211988.9, i0_err=827.0373, > chi2=0.447. > R1rho at 799.8 MHz, for offset=142.754 ppm and dispersion point 800.5, with 5 > time points. r2eff=3.240 r2eff_err=0.0253, i0=214417.4, i0_err=613.4378, > chi2=1.681. > R1rho at 799.8 MHz, for offset=142.754 ppm and dispersion point 1341.1, with > 5 time points. r2eff=5.084 r2eff_err=0.0352, i0=226358.7, i0_err=655.7670, > chi2=23.170. > R1rho at 799.8 MHz, for offset=179.768 ppm and dispersion point 1341.1, with > 5 time points. r2eff=2.208 r2eff_err=0.0178, i0=228620.6, i0_err=560.0873, > chi2=7.794. > R1rho at 799.8 MHz, for offset=241.459 ppm and dispersion point 1341.1, with > 5 time points. r2eff=1.711 r2eff_err=0.0155, i0=224087.5, i0_err=546.4217, > chi2=21.230. > > Fitting with minfx to: 52V @N > ----------------------------- > > min_algor='BFGS', c_code=False, constraints=False, chi2_jacobian?=True > ---------------------------------------------------------------------- > > R1rho at 799.8 MHz, for offset=118.078 ppm and dispersion point 431.0, with 4 > time points. r2eff=8.646 r2eff_err=0.0524, i0=202664.2, i0_err=1239.0827, > chi2=3.758. > R1rho at 799.8 MHz, for offset=118.078 ppm and dispersion point 651.2, with 5 > time points. r2eff=10.377 r2eff_err=0.0228, i0=206049.6, i0_err=178.1907, > chi2=27.291. > R1rho at 799.8 MHz, for offset=118.078 ppm and dispersion point 800.5, with 5 > time points. r2eff=10.506 r2eff_err=0.0345, i0=202586.3, i0_err=705.7630, > chi2=13.357. > R1rho at 799.8 MHz, for offset=118.078 ppm and dispersion point 984.0, with 5 > time points. r2eff=10.903 r2eff_err=0.0206, i0=203455.0, i0_err=186.0857, > chi2=33.632. > R1rho at 799.8 MHz, for offset=118.078 ppm and dispersion point 1341.1, with > 5 time points. r2eff=10.684 r2eff_err=0.0198, i0=218670.4, i0_err=165.0420, > chi2=35.818. > R1rho at 799.8 MHz, for offset=118.078 ppm and dispersion point 1648.5, with > 5 time points. r2eff=10.501 r2eff_err=0.0407, i0=206502.5, i0_err=321.3685, > chi2=7.356. > R1rho at 799.8 MHz, for offset=124.247 ppm and dispersion point 1341.1, with > 5 time points. r2eff=11.118 r2eff_err=0.0301, i0=216447.2, i0_err=248.9394, > chi2=15.587. > R1rho at 799.8 MHz, for offset=130.416 ppm and dispersion point 800.5, with 5 > time points. r2eff=7.866 r2eff_err=0.0280, i0=211869.7, i0_err=259.8845, > chi2=14.585. > R1rho at 799.8 MHz, for offset=130.416 ppm and dispersion point 1341.1, with > 5 time points. r2eff=9.259 r2eff_err=0.0108, i0=217703.2, i0_err=88.1514, > chi2=79.498. > R1rho at 799.8 MHz, for offset=130.416 ppm and dispersion point 1648.5, with > 5 time points. r2eff=9.565 r2eff_err=0.1630, i0=211988.9, i0_err=2054.6615, > chi2=0.447. > R1rho at 799.8 MHz, for offset=142.754 ppm and dispersion point 800.5, with 5 > time points. r2eff=3.240 r2eff_err=0.0485, i0=214417.4, i0_err=611.7573, > chi2=1.681. > R1rho at 799.8 MHz, for offset=142.754 ppm and dispersion point 1341.1, with > 5 time points. r2eff=5.084 r2eff_err=0.0124, i0=226358.7, i0_err=122.7341, > chi2=23.170. > R1rho at 799.8 MHz, for offset=179.768 ppm and dispersion point 1341.1, with > 5 time points. r2eff=2.208 r2eff_err=0.0086, i0=228620.6, i0_err=219.4208, > chi2=7.794. > R1rho at 799.8 MHz, for offset=241.459 ppm and dispersion point 1341.1, with > 5 time points. r2eff=1.711 r2eff_err=0.0101, i0=224087.5, i0_err=166.9081, > chi2=21.230. > > task #7822(https://gna.org/task/index.php?7822): Implement user function to > estimate R2eff and associated errors for exponential curve fitting. > > Modified: > trunk/specific_analyses/relax_disp/estimate_r2eff.py > trunk/test_suite/system_tests/relax_disp.py > > Modified: trunk/specific_analyses/relax_disp/estimate_r2eff.py > URL: > http://svn.gna.org/viewcvs/relax/trunk/specific_analyses/relax_disp/estimate_r2eff.py?rev=25379&r1=25378&r2=25379&view=diff > ============================================================================== > --- trunk/specific_analyses/relax_disp/estimate_r2eff.py (original) > +++ trunk/specific_analyses/relax_disp/estimate_r2eff.py Thu Aug 28 > 15:14:16 2014 > @@ -175,7 +175,7 @@ > print(print_string), > > > -def multifit_covar(J=None, epsrel=0.0, errors=None): > +def multifit_covar(J=None, epsrel=0.0, errors=None, use_weights=True): > """This is the implementation of the multifit covariance. > > This is inspired from GNU Scientific Library (GSL). > @@ -184,9 +184,15 @@ > > The parameter 'epsrel' is used to remove linear-dependent columns when J > is rank deficient. > > + The weighting matrix 'W', is a square symmetric matrix. For independent > measurements, this is a diagonal matrix. Larger values indicate greater > significance. It is formed by multiplying the supplied errors as 1./errors^2 > with an Identity matrix:: > + > + W = I.(1/errors^2) > + > + If 'use_weights' is set to 'False', the errors are set to 1.0. > + > The covariance matrix is given by:: > > - covar = (J^T J)^{-1} , > + covar = (J^T.W.J)^{-1} , > > and is computed by QR decomposition of J with column-pivoting. Any > columns of R which satisfy:: > > @@ -224,6 +230,8 @@ > @type epsrel: float > @keyword errors: The standard deviation of the measured intensity > values per time point. > @type errors: numpy array > + @keyword use_weights: If the supplied weights should be used. > + @type use_weights: bool > @return: The co-variance matrix > @rtype: square numpy array > """ > @@ -237,6 +245,10 @@ > # Now form the error matrix, with errors down the diagonal. > weights = 1. / errors**2 > > + if use_weights == False: > + weights[:] = 1.0 > + > + # Form weight matrix. > W = multiply(weights, eye_mat) > > # The covariance matrix (sometimes referred to as the > variance-covariance matrix), Qxx, is defined as: > @@ -344,7 +356,7 @@ > self.factor = factor > > > - def set_settings_minfx(self, scaling_matrix=None, min_algor='simplex', > c_code=True, constraints=False, func_tol=1e-25, grad_tol=None, > max_iterations=10000000): > + def set_settings_minfx(self, scaling_matrix=None, min_algor='simplex', > c_code=True, constraints=False, chi2_jacobian=False, func_tol=1e-25, > grad_tol=None, max_iterations=10000000): > """Setup options to minfx. > > @keyword scaling_matrix: The square and diagonal scaling matrix. > @@ -355,6 +367,8 @@ > @type c_code: bool > @keyword constraints: If constraints should be used. > @type constraints: bool > + @keyword chi2_jacobian: If the chi2 Jacobian should be used. > + @type chi2_jacobian: bool > @keyword func_tol: The function tolerance which, when > reached, terminates optimisation. Setting this to None turns of the check. > @type func_tol: None or float > @keyword grad_tol: The gradient tolerance which, when > reached, terminates optimisation. Setting this to None turns of the check. > @@ -366,6 +380,7 @@ > # Store variables. > self.scaling_matrix = scaling_matrix > self.c_code = c_code > + self.chi2_jacobian = chi2_jacobian > > # Scaling initialisation. > self.scaling_flag = False > @@ -561,7 +576,7 @@ > return 1. / self.errors * (self.func_exp(self.times, *params) - > self.values) > > > -def estimate_r2eff(method='minfx', min_algor='simplex', c_code=True, > constraints=False, spin_id=None, ftol=1e-15, xtol=1e-15, maxfev=10000000, > factor=100.0, verbosity=1): > +def estimate_r2eff(method='minfx', min_algor='simplex', c_code=True, > constraints=False, chi2_jacobian=False, spin_id=None, ftol=1e-15, xtol=1e-15, > maxfev=10000000, factor=100.0, verbosity=1): > """Estimate r2eff and errors by exponential curve fitting with > scipy.optimize.leastsq or minfx. > > THIS IS ONLY FOR TESTING. > @@ -583,10 +598,12 @@ > @type method: string > @keyword min_algor: The minimisation algorithm > @type min_algor: string > + @keyword c_code: If optimise with C code. > + @type c_code: bool > @keyword constraints: If constraints should be used. > @type constraints: bool > - @keyword c_code: If optimise with C code. > - @type c_code: bool > + @keyword chi2_jacobian: If the chi2 Jacobian should be used. > + @type chi2_jacobian: bool > @keyword spin_id: The spin identification string. > @type spin_id: str > @keyword ftol: The function tolerance for the relative > error desired in the sum of squares, parsed to leastsq. > @@ -661,7 +678,7 @@ > top += 2 > subsection(file=sys.stdout, text="Fitting with %s to: > %s"%(method, spin_string), prespace=top) > if method == 'minfx': > - subsection(file=sys.stdout, text="min_algor='%s', c_code=%s, > constraints=%s"%(min_algor, c_code, constraints), prespace=0) > + subsection(file=sys.stdout, text="min_algor='%s', c_code=%s, > constraints=%s, chi2_jacobian?=%s"%(min_algor, c_code, constraints, > chi2_jacobian), prespace=0) > > # Loop over each spectrometer frequency and dispersion point. > for exp_type, frq, offset, point, ei, mi, oi, di in > loop_exp_frq_offset_point(return_indices=True): > @@ -692,7 +709,7 @@ > > elif method == 'minfx': > # Set settings. > - E.set_settings_minfx(min_algor=min_algor, c_code=c_code, > constraints=constraints) > + E.set_settings_minfx(min_algor=min_algor, c_code=c_code, > chi2_jacobian=chi2_jacobian, constraints=constraints) > > # Acquire results. > results = minimise_minfx(E=E) > @@ -737,7 +754,7 @@ > point_info = "%s at %3.1f MHz, for offset=%3.3f ppm and > dispersion point %-5.1f, with %i time points." % (exp_type, frq/1E6, offset, > point, len(times)) > print_strings.append(point_info) > > - par_info = "r2eff=%3.3f r2eff_err=%3.3f, i0=%6.1f, > i0_err=%3.3f, chi2=%3.3f.\n" % ( r2eff, r2eff_err, i0, i0_err, chi2) > + par_info = "r2eff=%3.3f r2eff_err=%3.4f, i0=%6.1f, > i0_err=%3.4f, chi2=%3.3f.\n" % ( r2eff, r2eff_err, i0, i0_err, chi2) > print_strings.append(par_info) > > if E.verbosity >= 2: > @@ -912,14 +929,24 @@ > #jacobian_matrix_exp2 = E.jacobian_matrix_exp > #print jacobian_matrix_exp - jacobian_matrix_exp2 > else: > - # Call class, to store value. > - E.func_exp_grad(param_vector) > - jacobian_matrix_exp = E.jacobian_matrix_exp > - #E.func_exp_chi2_grad(param_vector) > - #jacobian_matrix_exp = E.jacobian_matrix_exp_chi2 > + if E.chi2_jacobian: > + # Call class, to store value. > + E.func_exp_chi2_grad(param_vector) > + jacobian_matrix_exp = E.jacobian_matrix_exp_chi2 > + else: > + # Call class, to store value. > + E.func_exp_grad(param_vector) > + jacobian_matrix_exp = E.jacobian_matrix_exp > + #E.func_exp_chi2_grad(param_vector) > + #jacobian_matrix_exp = E.jacobian_matrix_exp_chi2 > > # Get the co-variance > - pcov = multifit_covar(J=jacobian_matrix_exp, errors=E.errors) > + if E.chi2_jacobian: > + use_weights = False > + else: > + use_weights = True > + > + pcov = multifit_covar(J=jacobian_matrix_exp, errors=E.errors, > use_weights=use_weights) > > # To compute one standard deviation errors on the parameters, take the > square root of the diagonal covariance. > param_vector_error = sqrt(diag(pcov)) > > Modified: trunk/test_suite/system_tests/relax_disp.py > URL: > http://svn.gna.org/viewcvs/relax/trunk/test_suite/system_tests/relax_disp.py?rev=25379&r1=25378&r2=25379&view=diff > ============================================================================== > --- trunk/test_suite/system_tests/relax_disp.py (original) > +++ trunk/test_suite/system_tests/relax_disp.py Thu Aug 28 15:14:16 2014 > @@ -2946,12 +2946,13 @@ > > > # Now do it manually. > - estimate_r2eff(method='scipy.optimize.leastsq') > - estimate_r2eff(method='minfx', min_algor='simplex', c_code=True, > constraints=False) > - estimate_r2eff(method='minfx', min_algor='simplex', c_code=False, > constraints=False) > - estimate_r2eff(method='minfx', min_algor='BFGS', c_code=True, > constraints=False) > - estimate_r2eff(method='minfx', min_algor='BFGS', c_code=False, > constraints=False) > + #estimate_r2eff(method='scipy.optimize.leastsq') > + #estimate_r2eff(method='minfx', min_algor='simplex', c_code=True, > constraints=False) > + #estimate_r2eff(method='minfx', min_algor='simplex', c_code=False, > constraints=False) > + #estimate_r2eff(method='minfx', min_algor='BFGS', c_code=True, > constraints=False) > + #estimate_r2eff(method='minfx', min_algor='BFGS', c_code=False, > constraints=False) > estimate_r2eff(method='minfx', min_algor='Newton', c_code=True, > constraints=False) > + estimate_r2eff(method='minfx', min_algor='BFGS', c_code=False, > constraints=False, chi2_jacobian=True) > > > def test_exp_fit(self): > > > _______________________________________________ > relax (http://www.nmr-relax.com) > > This is the relax-commits mailing list > relax-comm...@gna.org > > To unsubscribe from this list, get a password > reminder, or change your subscription options, > visit the list information page at > https://mail.gna.org/listinfo/relax-commits _______________________________________________ relax (http://www.nmr-relax.com) This is the relax-devel mailing list relax-devel@gna.org To unsubscribe from this list, get a password reminder, or change your subscription options, visit the list information page at https://mail.gna.org/listinfo/relax-devel
