Could you create a simple relax script and log file in the test suite directories with the results of a super-massive Monte Carlo simulation run? I.e. a really insane number of Monte Carlo simulations and then a clear print out of both the R2eff and I0 errors. I really cannot follow the numbers in these tests. There are too many points where bugs could be hidden. It would be good to also have a simple relax script and log file for the covariance matrix estimate as well.
Cheers, Edward On 29 August 2014 11:59, Troels Emtekær Linnet <[email protected]> wrote: > You may want to look here: > > relax -s Relax_disp.test_estimate_r2eff_err_methods -d > > 2014-08-29 11:57 GMT+02:00 Troels Emtekær Linnet <[email protected]>: >> Hi Edward. >> >> There is something totally wrong with the C, Jacobian. >> Errors are estimated to: >> >> 37.619 17.290 25.616 16.036 16.164 32.826 22.920 21.462 7.777 145.309 >> 36.884 9.116 6.199 7.018 sum= 402.235 >> >> Which is much different to: >> 0.041 0.040 0.040 0.054 0.041 0.044 0.042 0.037 0.034 0.043 0.013 >> 0.018 0.007 0.010 sum= 0.462 >> >> You can see how the error estimation develops in: >> verify_estimate_r2eff_err_compare_mc >> >> You will see, that just 50 monte carlo simulations is better than estimating. >> >> Best >> Troels >> >> >> 2014-08-29 11:51 GMT+02:00 Edward d'Auvergne <[email protected]>: >>> 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
