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 <tlin...@nmr-relax.com>:
> 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 <edw...@nmr-relax.com>:
>> 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 <tlin...@nmr-relax.com> 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 <edw...@nmr-relax.com>:
>>>> 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, <tlin...@nmr-relax.com> 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
>>>>> 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
