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 <edward@xxxxxxxxxxxxx>:
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, <tlinnet@xxxxxxxxxxxxx> 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-commits@xxxxxxx
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@xxxxxxx
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
Re: r25411 - in /trunk: specific_analyses/relax_disp/ test_suite/system_tests/ user_functions/