Hi, I saw the results from that 'hidden' system test and was wondering what was happening? The Jacobian of the chi-squared function should remove the factor of 2, as it has a factor of minus two. But it also includes the difference between the measured and back-calculated peak intensities divided by the variance as well. So why does this Jacobian, which is much closer to the 2000 MC simulations, not work? I cannot understand this as it is totally illogical. If your error estimate is closer to the real thing, then you should get closer to the real optimisation results.
Do you have a log file somewhere which contains the results from the 2000 MC simulations? It might be worth creating a file which compares this, or even more simulations, 100,000 for example, to the covariance technique. Once the error estimate technique is functional and debugged, then we can work out why the models are optimisating differently. These two problems need to be separated and solved independently, otherwise you can encounter the common yet fatal coding problem of two opposing bugs partially cancelling out their effects. Regards, Edward On 29 August 2014 11:01, Troels Emtekær Linnet <[email protected]> wrote: > Hi Edward. > > Would it be possible to have both? > > The exponential Jacobian, and the chi2 Jacobian. > > My tests last night showed something weird. > > Using the chi2 Jacobian, the errors come closer to the ones reported > my MC calculations. > The direct jacobian would have double error on R2eff. > > But when fitting for R1rho models, using the errors from the direct > jacobian, was much better in agreement with > MC error fitting. > > The parameters from chi2 Jacobian, was worse. > > See verify_r1rho_kjaergaard_missing_r1() in systemtest for comparison. > > Look at the 'kex' parameter! > > # Compare values. > if spin_id == ':52@N': > if param == 'r1': > if model == MODEL_NOREX: > if r2eff_estimate == 'direct': > self.assertAlmostEqual(value, 1.46138805) > elif r2eff_estimate == 'MC2000': > self.assertAlmostEqual(value, 1.46328102) > elif r2eff_estimate == 'chi2': > self.assertAlmostEqual(value, 1.43820629) > elif model == MODEL_DPL94: > if r2eff_estimate == 'direct': > self.assertAlmostEqual(value, 1.44845742) > elif r2eff_estimate == 'MC2000': > self.assertAlmostEqual(value, 1.45019848) > elif r2eff_estimate == 'chi2': > self.assertAlmostEqual(value, 1.44666512) > elif model == MODEL_TP02: > if r2eff_estimate == 'direct': > self.assertAlmostEqual(value, 1.54354392) > elif r2eff_estimate == 'MC2000': > self.assertAlmostEqual(value, 1.54352369) > elif r2eff_estimate == 'chi2': > self.assertAlmostEqual(value, 1.55964020) > elif model == MODEL_TAP03: > if r2eff_estimate == 'direct': > self.assertAlmostEqual(value, 1.54356410) > elif r2eff_estimate == 'MC2000': > self.assertAlmostEqual(value, 1.54354367) > elif r2eff_estimate == 'chi2': > self.assertAlmostEqual(value, 1.55967157) > elif model == MODEL_MP05: > if r2eff_estimate == 'direct': > self.assertAlmostEqual(value, 1.54356416) > elif r2eff_estimate == 'MC2000': > self.assertAlmostEqual(value, 1.54354372) > elif r2eff_estimate == 'chi2': > self.assertAlmostEqual(value, 1.55967163) > elif model == MODEL_NS_R1RHO_2SITE: > if r2eff_estimate == 'direct': > self.assertAlmostEqual(value, 1.41359221, 5) > elif r2eff_estimate == 'MC2000': > self.assertAlmostEqual(value, 1.41321968, 5) > elif r2eff_estimate == 'chi2': > self.assertAlmostEqual(value, 1.36303129, 5) > > elif param == 'r2': > if model == MODEL_NOREX: > if r2eff_estimate == 'direct': > self.assertAlmostEqual(value, 11.48392439) > elif r2eff_estimate == 'MC2000': > self.assertAlmostEqual(value, 11.48040934) > elif r2eff_estimate == 'chi2': > self.assertAlmostEqual(value, 11.47224488) > elif model == MODEL_DPL94: > if r2eff_estimate == 'direct': > self.assertAlmostEqual(value, 10.15688372, 6) > elif r2eff_estimate == 'MC2000': > self.assertAlmostEqual(value, 10.16304887, 6) > elif r2eff_estimate == 'chi2': > self.assertAlmostEqual(value, 9.20037797, 6) > elif model == MODEL_TP02: > if r2eff_estimate == 'direct': > self.assertAlmostEqual(value, 9.72654896, 6) > elif r2eff_estimate == 'MC2000': > self.assertAlmostEqual(value, 9.72772726, 6) > elif r2eff_estimate == 'chi2': > self.assertAlmostEqual(value, 9.53948340, 6) > elif model == MODEL_TAP03: > if r2eff_estimate == 'direct': > self.assertAlmostEqual(value, 9.72641887, 6) > elif r2eff_estimate == 'MC2000': > self.assertAlmostEqual(value, 9.72759374, 6) > elif r2eff_estimate == 'chi2': > self.assertAlmostEqual(value, 9.53926913, 6) > elif model == MODEL_MP05: > if r2eff_estimate == 'direct': > self.assertAlmostEqual(value, 9.72641723, 6) > elif r2eff_estimate == 'MC2000': > self.assertAlmostEqual(value, 9.72759220, 6) > elif r2eff_estimate == 'chi2': > self.assertAlmostEqual(value, 9.53926778, 6) > elif model == MODEL_NS_R1RHO_2SITE: > if r2eff_estimate == 'direct': > self.assertAlmostEqual(value, 9.34531535, 5) > elif r2eff_estimate == 'MC2000': > self.assertAlmostEqual(value, 9.34602793, 5) > elif r2eff_estimate == 'chi2': > self.assertAlmostEqual(value, 9.17631409, 5) > > # For all other parameters. > else: > # Get the value. > value = getattr(cur_spin, param) > > # Print value. > print("%-10s %-6s %-6s %3.8f" % ("Parameter:", param, "Value:", value)) > > # Compare values. > if spin_id == ':52@N': > if param == 'phi_ex': > if model == MODEL_DPL94: > if r2eff_estimate == 'direct': > self.assertAlmostEqual(value, 0.07599563) > elif r2eff_estimate == 'MC2000': > self.assertAlmostEqual(value, 0.07561937) > elif r2eff_estimate == 'chi2': > self.assertAlmostEqual(value, 0.12946061) > > elif param == 'pA': > if model == MODEL_TP02: > if r2eff_estimate == 'direct': > self.assertAlmostEqual(value, 0.88827040) > elif r2eff_estimate == 'MC2000': > self.assertAlmostEqual(value, 0.88807487) > elif r2eff_estimate == 'chi2': > self.assertAlmostEqual(value, 0.87746233) > elif model == MODEL_TAP03: > if r2eff_estimate == 'direct': > self.assertAlmostEqual(value, 0.88828922) > elif r2eff_estimate == 'MC2000': > self.assertAlmostEqual(value, 0.88809318) > elif r2eff_estimate == 'chi2': > self.assertAlmostEqual(value, 0.87747558) > elif model == MODEL_MP05: > if r2eff_estimate == 'direct': > self.assertAlmostEqual(value, 0.88828924) > elif r2eff_estimate == 'MC2000': > self.assertAlmostEqual(value, 0.88809321) > elif r2eff_estimate == 'chi2': > self.assertAlmostEqual(value, 0.87747562) > elif model == MODEL_NS_R1RHO_2SITE: > if r2eff_estimate == 'direct': > self.assertAlmostEqual(value, 0.94504369, 6) > elif r2eff_estimate == 'MC2000': > self.assertAlmostEqual(value, 0.94496541, 6) > elif r2eff_estimate == 'chi2': > self.assertAlmostEqual(value, 0.92084707, 6) > > elif param == 'dw': > if model == MODEL_TP02: > if r2eff_estimate == 'direct': > self.assertAlmostEqual(value, 1.08875840, 6) > elif r2eff_estimate == 'MC2000': > self.assertAlmostEqual(value, 1.08765638, 6) > elif r2eff_estimate == 'chi2': > self.assertAlmostEqual(value, 1.09753230, 6) > elif model == MODEL_TAP03: > if r2eff_estimate == 'direct': > self.assertAlmostEqual(value, 1.08837238, 6) > elif r2eff_estimate == 'MC2000': > self.assertAlmostEqual(value, 1.08726698, 6) > elif r2eff_estimate == 'chi2': > self.assertAlmostEqual(value, 1.09708821, 6) > elif model == MODEL_MP05: > if r2eff_estimate == 'direct': > self.assertAlmostEqual(value, 1.08837241, 6) > elif r2eff_estimate == 'MC2000': > self.assertAlmostEqual(value, 1.08726706, 6) > elif r2eff_estimate == 'chi2': > self.assertAlmostEqual(value, 1.09708832, 6) > elif model == MODEL_NS_R1RHO_2SITE: > if r2eff_estimate == 'direct': > self.assertAlmostEqual(value, 1.56001812, 5) > elif r2eff_estimate == 'MC2000': > self.assertAlmostEqual(value, 1.55833321, 5) > elif r2eff_estimate == 'chi2': > self.assertAlmostEqual(value, 1.36406712, 5) > > elif param == 'kex': > if model == MODEL_DPL94: > if r2eff_estimate == 'direct': > self.assertAlmostEqual(value, 4460.43711569, 2) > elif r2eff_estimate == 'MC2000': > self.assertAlmostEqual(value, 4419.03917195, 2) > elif r2eff_estimate == 'chi2': > self.assertAlmostEqual(value, 6790.22736344, 2) > elif model == MODEL_TP02: > if r2eff_estimate == 'direct': > self.assertAlmostEqual(value, 4921.28602757, 3) > elif r2eff_estimate == 'MC2000': > self.assertAlmostEqual(value, 4904.70144883, 3) > elif r2eff_estimate == 'chi2': > self.assertAlmostEqual(value, 5146.20306591, 3) > elif model == MODEL_TAP03: > if r2eff_estimate == 'direct': > self.assertAlmostEqual(value, 4926.42963491, 3) > elif r2eff_estimate == 'MC2000': > self.assertAlmostEqual(value, 4909.86877150, 3) > elif r2eff_estimate == 'chi2': > self.assertAlmostEqual(value, 5152.51105814, 3) > elif model == MODEL_MP05: > if r2eff_estimate == 'direct': > self.assertAlmostEqual(value, 4926.44236315, 3) > elif r2eff_estimate == 'MC2000': > self.assertAlmostEqual(value, 4909.88110195, 3) > elif r2eff_estimate == 'chi2': > self.assertAlmostEqual(value, 5152.52097111, 3) > elif model == MODEL_NS_R1RHO_2SITE: > if r2eff_estimate == 'direct': > self.assertAlmostEqual(value, 5628.66061488, 2) > elif r2eff_estimate == 'MC2000': > self.assertAlmostEqual(value, 5610.20221435, 2) > elif r2eff_estimate == 'chi2': > self.assertAlmostEqual(value, 5643.34067090, 2) > > elif param == 'chi2': > if model == MODEL_NOREX: > if r2eff_estimate == 'direct': > self.assertAlmostEqual(value, 848.42016907, 5) > elif r2eff_estimate == 'MC2000': > self.assertAlmostEqual(value, 3363.95829122, 5) > elif r2eff_estimate == 'chi2': > self.assertAlmostEqual(value, 5976.49946726, 5) > elif model == MODEL_DPL94: > if r2eff_estimate == 'direct': > self.assertAlmostEqual(value, 179.47041241) > elif r2eff_estimate == 'MC2000': > self.assertAlmostEqual(value, 710.24767560) > elif r2eff_estimate == 'chi2': > self.assertAlmostEqual(value, 612.72616697, 5) > elif model == MODEL_TP02: > if r2eff_estimate == 'direct': > self.assertAlmostEqual(value, 29.33882530, 6) > elif r2eff_estimate == 'MC2000': > self.assertAlmostEqual(value, 114.47142772, 6) > elif r2eff_estimate == 'chi2': > self.assertAlmostEqual(value, 250.50838162, 5) > elif model == MODEL_TAP03: > if r2eff_estimate == 'direct': > self.assertAlmostEqual(value, 29.29050673, 6) > elif r2eff_estimate == 'MC2000': > self.assertAlmostEqual(value, 114.27987534) > elif r2eff_estimate == 'chi2': > self.assertAlmostEqual(value, 250.04050719, 5) > elif model == MODEL_MP05: > if r2eff_estimate == 'direct': > self.assertAlmostEqual(value, 29.29054301, 6) > elif r2eff_estimate == 'MC2000': > self.assertAlmostEqual(value, 114.28002272) > elif r2eff_estimate == 'chi2': > self.assertAlmostEqual(value, 250.04077478, 5) > elif model == MODEL_NS_R1RHO_2SITE: > if r2eff_estimate == 'direct': > self.assertAlmostEqual(value, 34.44010543, 6) > elif r2eff_estimate == 'MC2000': > self.assertAlmostEqual(value, 134.14368365) > elif r2eff_estimate == 'chi2': > self.assertAlmostEqual(value, 278.55121388, 5) > > 2014-08-29 9:49 GMT+02:00 Edward d'Auvergne <[email protected]>: >> Hi Troels, >> >> I've now converted the target_functions.relax_fit.jacobian() function >> to be the Jacobian of the chi-squared function rather than the >> Jacobian of the exponential function. This should match your >> specific_analyses.relax_disp.estimate_r2eff.func_exp_chi2_grad() >> function. I mixed up the two because the Levenberg-Marquardt >> algorithm in minfx requires the Jacobian of the exponential, and it's >> been 8 years since I last derived and implemented a Jacobian. >> >> Regards, >> >> Edward >> >> >> >> On 28 August 2014 21:43, <[email protected]> wrote: >>> Author: tlinnet >>> Date: Thu Aug 28 21:43:13 2014 >>> New Revision: 25411 >>> >>> URL: http://svn.gna.org/viewcvs/relax?rev=25411&view=rev >>> Log: >>> Reverted the logic, that the chi2 Jacobian should be used. >>> >>> Instead, the direct Jacobian exponential is used instead. >>> >>> When fitting with the estimated errors from the Direct Jacobian, the >>> results are MUCH better, and comparable >>> to 2000 Monte-Carlo simulations. >>> >>> 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 >>> trunk/user_functions/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=25411&r1=25410&r2=25411&view=diff >>> ============================================================================== >>> --- trunk/specific_analyses/relax_disp/estimate_r2eff.py (original) >>> +++ trunk/specific_analyses/relax_disp/estimate_r2eff.py Thu Aug 28 >>> 21:43:13 2014 >>> @@ -90,7 +90,7 @@ >>> return jacobian_matrix_exp_chi2 >>> >>> >>> -def estimate_r2eff_err(chi2_jacobian=True, spin_id=None, epsrel=0.0, >>> verbosity=1): >>> +def estimate_r2eff_err(chi2_jacobian=False, spin_id=None, epsrel=0.0, >>> verbosity=1): >>> """This will estimate the R2eff and i0 errors from the covariance >>> matrix Qxx. Qxx is calculated from the Jacobian matrix and the optimised >>> parameters. >>> >>> @keyword chi2_jacobian: If the Jacobian derived from the chi2 >>> function, should be used instead of the Jacobian from the exponential >>> function. >>> >>> 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=25411&r1=25410&r2=25411&view=diff >>> ============================================================================== >>> --- trunk/test_suite/system_tests/relax_disp.py (original) >>> +++ trunk/test_suite/system_tests/relax_disp.py Thu Aug 28 21:43:13 2014 >>> @@ -2744,13 +2744,13 @@ >>> self.interpreter.minimise.execute(min_algor='Newton', >>> constraints=False, verbosity=1) >>> >>> # Estimate R2eff errors. >>> - self.interpreter.relax_disp.r2eff_err_estimate(chi2_jacobian=False) >>> + self.interpreter.relax_disp.r2eff_err_estimate(chi2_jacobian=True) >>> >>> # Run the analysis. >>> relax_disp.Relax_disp(pipe_name=ds.pipe_name, >>> pipe_bundle=ds.pipe_bundle, results_dir=result_dir_name, models=MODELS, >>> grid_inc=GRID_INC, mc_sim_num=MC_NUM, modsel=MODSEL) >>> >>> # Verify the data. >>> - self.verify_r1rho_kjaergaard_missing_r1(models=MODELS, >>> result_dir_name=result_dir_name, r2eff_estimate='direct') >>> + self.verify_r1rho_kjaergaard_missing_r1(models=MODELS, >>> result_dir_name=result_dir_name, r2eff_estimate='chi2') >>> >>> >>> def test_estimate_r2eff_err_auto(self): >>> @@ -2849,7 +2849,7 @@ >>> relax_disp.Relax_disp(pipe_name=pipe_name, >>> pipe_bundle=pipe_bundle, results_dir=result_dir_name, models=MODELS, >>> grid_inc=GRID_INC, mc_sim_num=MC_NUM, exp_mc_sim_num=EXP_MC_NUM, >>> modsel=MODSEL, r1_fit=r1_fit) >>> >>> # Verify the data. >>> - self.verify_r1rho_kjaergaard_missing_r1(models=MODELS, >>> result_dir_name=result_dir_name, r2eff_estimate='chi2') >>> + self.verify_r1rho_kjaergaard_missing_r1(models=MODELS, >>> result_dir_name=result_dir_name, r2eff_estimate='direct') >>> >>> >>> def test_estimate_r2eff_err_methods(self): >>> >>> Modified: trunk/user_functions/relax_disp.py >>> URL: >>> http://svn.gna.org/viewcvs/relax/trunk/user_functions/relax_disp.py?rev=25411&r1=25410&r2=25411&view=diff >>> ============================================================================== >>> --- trunk/user_functions/relax_disp.py (original) >>> +++ trunk/user_functions/relax_disp.py Thu Aug 28 21:43:13 2014 >>> @@ -636,7 +636,7 @@ >>> uf.title_short = "Estimate R2eff errors." >>> uf.add_keyarg( >>> name = "chi2_jacobian", >>> - default = True, >>> + default = False, >>> py_type = "bool", >>> desc_short = "use of chi2 Jacobian", >>> desc = "If the Jacobian derived from the chi2 function, should be used >>> instead of the Jacobian from the exponential function." >>> >>> >>> _______________________________________________ >>> relax (http://www.nmr-relax.com) >>> >>> This is the relax-commits mailing list >>> [email protected] >>> >>> 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 >> [email protected] >> >> 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 [email protected] 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
