Hi, No, that's not correct. Try performing the maths yourself and try to derive the chi-squared partial derivative. You will see that it's a little different.
Regards, Edward On 28 August 2014 15:38, Troels Emtekær Linnet <tlin...@nmr-relax.com> wrote: > Hi Edward. > > I is just target_functions/c_chi2.c where you dont sum the elements, but > return as array. > > Best > Troels > > 2014-08-28 15:30 GMT+02:00 Edward d'Auvergne <edw...@nmr-relax.com>: >> 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 _______________________________________________ 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
