Hi Edward.
The new Jacobian, does not do the Job:
-2, is the implemented C code chi2 Jacobian.
-1, is:
-------
# Make partial derivative, with respect to r2eff.
d_chi2_d_r2eff = 2.0 * i0 * times * ( -i0 * exp( -r2eff * times) +
values) * exp( -r2eff * times ) / errors**2
# Make partial derivative, with respect to i0.
d_chi2_d_i0 = - 2.0 * ( -i0 * exp( -r2eff * times) + values) *
exp( -r2eff * times) / errors**2
---------
Results from:
Relax_disp.verify_estimate_r2eff_err_compare_mc
-2 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
-1 0.052 0.023 0.034 0.021 0.020 0.041 0.030 0.028 0.011 0.163
0.048 0.012 0.009 0.010 sum= 0.502
0 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
0.000 0.000 0.000 0.000 sum= 0.000
10 0.041 0.052 0.040 0.042 0.040 0.037 0.027 0.017 0.027 0.043
0.008 0.015 0.010 0.006 sum= 0.405
20 0.029 0.036 0.037 0.040 0.040 0.028 0.050 0.035 0.035 0.042
0.017 0.018 0.010 0.007 sum= 0.424
30 0.036 0.038 0.036 0.054 0.042 0.036 0.049 0.034 0.034 0.038
0.014 0.018 0.010 0.008 sum= 0.447
40 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
2014-08-29 11:01 GMT+02:00 Troels Emtekær Linnet <[email protected]>:
> 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
