Hi, That is not good :S With my experience with Scipy optimisation, you could probably guess that I would blame a bug in their covariance matrix code ;) In any case, this is clearly a bug and it needs to be solved to have useful R2eff errors from the covariance technique. Maybe the new Jacobian code could be used to create an independent covariance matrix to compare to:
http://thread.gmane.org/gmane.science.nmr.relax.scm/22980/focus=6828 Though the target_functions.relax_fit.jacobian() code may still need debugging, it could be a useful aid for debugging the Scipy covariance matrix. Regards, Edward On 25 August 2014 16:53, Troels Emtekær Linnet <tlin...@nmr-relax.com> wrote: > Hi Edward. > > The error for r2eff seems twice as high, while i0_err are almost identical. > > Best > Troels > > R1rho at 799.8 MHz, for offset=118.078 ppm and dispersion point 431.0. > r2eff=8.646/8.646 r2eff_err=0.0348/0.0692 i0=202664.191/202664.191 > i0_err=699.6443/712.4201 > > R1rho at 799.8 MHz, for offset=118.078 ppm and dispersion point 651.2. > r2eff=10.377/10.377 r2eff_err=0.0403/0.0810 i0=206049.558/206049.558 > i0_err=776.4215/782.1833 > > R1rho at 799.8 MHz, for offset=118.078 ppm and dispersion point 800.5. > r2eff=10.506/10.506 r2eff_err=0.0440/0.0853 i0=202586.332/202586.332 > i0_err=763.9678/758.7052 > > R1rho at 799.8 MHz, for offset=118.078 ppm and dispersion point 984.0. > r2eff=10.903/10.903 r2eff_err=0.0476/0.0922 i0=203455.021/203455.021 > i0_err=837.8779/828.7280 > > R1rho at 799.8 MHz, for offset=118.078 ppm and dispersion point 1341.1. > r2eff=10.684/10.684 r2eff_err=0.0446/0.0853 i0=218670.412/218670.412 > i0_err=850.0210/830.9558 > > R1rho at 799.8 MHz, for offset=118.078 ppm and dispersion point 1648.5. > r2eff=10.501/10.501 r2eff_err=0.0371/0.0742 i0=206502.512/206502.512 > i0_err=794.0523/772.9843 > > R1rho at 799.8 MHz, for offset=124.247 ppm and dispersion point 1341.1. > r2eff=11.118/11.118 r2eff_err=0.0413/0.0827 i0=216447.241/216447.241 > i0_err=784.6562/788.0384 > > R1rho at 799.8 MHz, for offset=130.416 ppm and dispersion point 800.5. > r2eff=7.866/7.866 r2eff_err=0.0347/0.0695 i0=211869.715/211869.715 > i0_err=749.2776/763.6930 > > R1rho at 799.8 MHz, for offset=130.416 ppm and dispersion point 1341.1. > r2eff=9.259/9.259 r2eff_err=0.0331/0.0661 i0=217703.151/217703.151 > i0_err=682.2137/685.5838 > > R1rho at 799.8 MHz, for offset=130.416 ppm and dispersion point 1648.5. > r2eff=9.565/9.565 r2eff_err=0.0373/0.0745 i0=211988.939/211988.939 > i0_err=839.0313/827.0373 > > R1rho at 799.8 MHz, for offset=142.754 ppm and dispersion point 800.5. > r2eff=3.240/3.240 r2eff_err=0.0127/0.0253 i0=214417.382/214417.382 > i0_err=595.8865/613.4378 > > R1rho at 799.8 MHz, for offset=142.754 ppm and dispersion point 1341.1. > r2eff=5.084/5.084 r2eff_err=0.0177/0.0352 i0=226358.691/226358.691 > i0_err=660.5314/655.7670 > > R1rho at 799.8 MHz, for offset=179.768 ppm and dispersion point 1341.1. > r2eff=2.208/2.208 r2eff_err=0.0091/0.0178 i0=228620.553/228620.553 > i0_err=564.8353/560.0873 > > R1rho at 799.8 MHz, for offset=241.459 ppm and dispersion point 1341.1. > r2eff=1.711/1.711 r2eff_err=0.0077/0.0155 i0=224087.486/224087.486 > i0_err=539.4300/546.4217 > > 2014-08-25 16:35 GMT+02:00 Edward d'Auvergne <edw...@nmr-relax.com>: >> Hi, >> >> The errors from the covariance matrix estimate, via >> scipy.optimize.leastsq() should be similar or even identical to the >> Monte Carlo simulation errors, not double. Is this still the case? >> >> Regards, >> >> Edward >> >> On 25 August 2014 01:32, <tlin...@nmr-relax.com> wrote: >>> Author: tlinnet >>> Date: Mon Aug 25 01:32:13 2014 >>> New Revision: 25238 >>> >>> URL: http://svn.gna.org/viewcvs/relax?rev=25238&view=rev >>> Log: >>> Modified systemtest Relax_disp.test_estimate_r2eff. >>> >>> This is to compare against errors simulated with 2000 MC. >>> >>> The paramaters are comparable, but not equal. >>> Mostly, it seems that the errors from scipy.optimize.leastsq, are twice as >>> high than the Monte Carlo simulations. >>> This affect model fitting, and the calculated chi2 value. >>> >>> Left column is 2000 Monte Carlo, right column is scipy.optimize.leastsq. >>> >>> Optimised parameters for spin: 52V @N >>> Model: No Rex >>> Parameter: r1 Value: 1.46138806 - 1.46328102 >>> Parameter: r2 Value: 11.48392438 - 11.48040934 >>> Parameter: chi2 Value: 848.42015672 - 3363.95829122 >>> >>> Model: DPL94 >>> Parameter: r1 Value: 1.44845743 - 1.45019848 >>> Parameter: r2 Value: 10.15688373 - 10.16304892 >>> Parameter: phi_ex Value: 0.07599563 - 0.07561937 >>> Parameter: kex Value: 4460.43707304 - 4419.03906628 >>> Parameter: chi2 Value: 179.47041255 - 710.24767560 >>> >>> Model: TP02 >>> Parameter: r1 Value: 1.54354392 - 1.54352369 >>> Parameter: r2 Value: 9.72654895 - 9.72772727 >>> Parameter: pA Value: 0.88827039 - 0.88807488 >>> Parameter: dw Value: 1.08875836 - 1.08765645 >>> Parameter: kex Value: 4921.28597928 - 4904.70134941 >>> Parameter: chi2 Value: 29.33882481 - 114.47142772 >>> >>> Model: TAP03 >>> Parameter: r1 Value: 1.54356410 - 1.54354368 >>> Parameter: r2 Value: 9.72641885 - 9.72759371 >>> Parameter: pA Value: 0.88828925 - 0.88809317 >>> Parameter: dw Value: 1.08837248 - 1.08726695 >>> Parameter: kex Value: 4926.42974479 - 4909.86896567 >>> Parameter: chi2 Value: 29.29050624 - 114.27987534 >>> >>> Model: MP05 >>> Parameter: r1 Value: 1.54356415 - 1.54354372 >>> Parameter: r2 Value: 9.72641730 - 9.72759220 >>> Parameter: pA Value: 0.88828927 - 0.88809322 >>> Parameter: dw Value: 1.08837250 - 1.08726707 >>> Parameter: kex Value: 4926.44228958 - 4909.88128236 >>> Parameter: chi2 Value: 29.29054252 - 114.28002272 >>> >>> Model: NS R1rho 2-site >>> Parameter: r1 Value: 1.41359226 - 1.41321968 >>> Parameter: r2 Value: 9.34531364 - 9.34602793 >>> Parameter: pA Value: 0.94504369 - 0.94496541 >>> Parameter: dw Value: 1.56001843 - 1.55833321 >>> Parameter: kex Value: 5628.66529504 - 5610.20221435 >>> Parameter: chi2 Value: 34.44010458 - 134.14368365 >>> >>> task #7822(https://gna.org/task/index.php?7822): Implement user function to >>> estimate R2eff and associated errors for exponential curve fitting. >>> >>> Modified: >>> trunk/test_suite/system_tests/relax_disp.py >>> >>> 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=25238&r1=25237&r2=25238&view=diff >>> ============================================================================== >>> --- trunk/test_suite/system_tests/relax_disp.py (original) >>> +++ trunk/test_suite/system_tests/relax_disp.py Mon Aug 25 01:32:13 2014 >>> @@ -2679,7 +2679,7 @@ >>> self.setup_r1rho_kjaergaard(cluster_ids=cluster_ids, read_R1=False) >>> >>> # The dispersion models. >>> - MODELS = [MODEL_NOREX, MODEL_DPL94] >>> + MODELS = [MODEL_NOREX, MODEL_DPL94, MODEL_TP02, MODEL_TAP03, >>> MODEL_MP05, MODEL_NS_R1RHO_2SITE] >>> >>> # The grid search size (the number of increments per dimension). >>> GRID_INC = None >>> @@ -2723,6 +2723,9 @@ >>> >>> # 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, do_assert=False) >>> >>> >>> def test_exp_fit(self): >>> @@ -7343,7 +7346,7 @@ >>> w_eff_file.close() >>> >>> >>> - def verify_r1rho_kjaergaard_missing_r1(self, models=None, >>> result_dir_name=None): >>> + def verify_r1rho_kjaergaard_missing_r1(self, models=None, >>> result_dir_name=None, do_assert=True): >>> """Verification of test_r1rho_kjaergaard_missing_r1.""" >>> >>> # Check the kex value of residue 52 >>> @@ -7382,35 +7385,36 @@ >>> # Print value. >>> print("%-10s %-6s %-6s %3.8f" % ("Parameter:", >>> param, "Value:", value)) >>> >>> - # Compare values. >>> - if spin_id == ':52@N': >>> - if param == 'r1': >>> - if model == MODEL_NOREX: >>> - self.assertAlmostEqual(value, >>> 1.46328102) >>> - elif model == MODEL_DPL94: >>> - self.assertAlmostEqual(value, >>> 1.45019848) >>> - elif model == MODEL_TP02: >>> - self.assertAlmostEqual(value, >>> 1.54352369) >>> - elif model == MODEL_TAP03: >>> - self.assertAlmostEqual(value, >>> 1.54354367) >>> - elif model == MODEL_MP05: >>> - self.assertAlmostEqual(value, >>> 1.54354372) >>> - elif model == MODEL_NS_R1RHO_2SITE: >>> - self.assertAlmostEqual(value, >>> 1.41321968, 6) >>> - >>> - elif param == 'r2': >>> - if model == MODEL_NOREX: >>> - self.assertAlmostEqual(value, >>> 11.48040934) >>> - elif model == MODEL_DPL94: >>> - self.assertAlmostEqual(value, >>> 10.16304887, 6) >>> - elif model == MODEL_TP02: >>> - self.assertAlmostEqual(value, >>> 9.72772726, 6) >>> - elif model == MODEL_TAP03: >>> - self.assertAlmostEqual(value, >>> 9.72759374, 6) >>> - elif model == MODEL_MP05: >>> - self.assertAlmostEqual(value, >>> 9.72759220, 6) >>> - elif model == MODEL_NS_R1RHO_2SITE: >>> - self.assertAlmostEqual(value, >>> 9.34602793, 5) >>> + if do_assert: >>> + # Compare values. >>> + if spin_id == ':52@N': >>> + if param == 'r1': >>> + if model == MODEL_NOREX: >>> + self.assertAlmostEqual(value, >>> 1.46328102) >>> + elif model == MODEL_DPL94: >>> + self.assertAlmostEqual(value, >>> 1.45019848) >>> + elif model == MODEL_TP02: >>> + self.assertAlmostEqual(value, >>> 1.54352369) >>> + elif model == MODEL_TAP03: >>> + self.assertAlmostEqual(value, >>> 1.54354367) >>> + elif model == MODEL_MP05: >>> + self.assertAlmostEqual(value, >>> 1.54354372) >>> + elif model == MODEL_NS_R1RHO_2SITE: >>> + self.assertAlmostEqual(value, >>> 1.41321968, 6) >>> + >>> + elif param == 'r2': >>> + if model == MODEL_NOREX: >>> + self.assertAlmostEqual(value, >>> 11.48040934) >>> + elif model == MODEL_DPL94: >>> + self.assertAlmostEqual(value, >>> 10.16304887, 6) >>> + elif model == MODEL_TP02: >>> + self.assertAlmostEqual(value, >>> 9.72772726, 6) >>> + elif model == MODEL_TAP03: >>> + self.assertAlmostEqual(value, >>> 9.72759374, 6) >>> + elif model == MODEL_MP05: >>> + self.assertAlmostEqual(value, >>> 9.72759220, 6) >>> + elif model == MODEL_NS_R1RHO_2SITE: >>> + self.assertAlmostEqual(value, >>> 9.34602793, 5) >>> >>> # For all other parameters. >>> else: >>> @@ -7420,57 +7424,58 @@ >>> # 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: >>> - self.assertAlmostEqual(value, >>> 0.07561937) >>> - >>> - elif param == 'pA': >>> - if model == MODEL_TP02: >>> - self.assertAlmostEqual(value, >>> 0.88807487) >>> - elif model == MODEL_TAP03: >>> - self.assertAlmostEqual(value, >>> 0.88809318) >>> - elif model == MODEL_MP05: >>> - self.assertAlmostEqual(value, >>> 0.88809321) >>> - elif model == MODEL_NS_R1RHO_2SITE: >>> - self.assertAlmostEqual(value, >>> 0.94496541, 6) >>> - >>> - elif param == 'dw': >>> - if model == MODEL_TP02: >>> - self.assertAlmostEqual(value, >>> 1.08765638, 6) >>> - elif model == MODEL_TAP03: >>> - self.assertAlmostEqual(value, >>> 1.08726698, 6) >>> - elif model == MODEL_MP05: >>> - self.assertAlmostEqual(value, >>> 1.08726706, 6) >>> - elif model == MODEL_NS_R1RHO_2SITE: >>> - self.assertAlmostEqual(value, >>> 1.55833321, 5) >>> - >>> - elif param == 'kex': >>> - if model == MODEL_DPL94: >>> - self.assertAlmostEqual(value, >>> 4419.03917195, 2) >>> - elif model == MODEL_TP02: >>> - self.assertAlmostEqual(value, >>> 4904.70144883, 3) >>> - elif model == MODEL_TAP03: >>> - self.assertAlmostEqual(value, >>> 4909.86877150, 3) >>> - elif model == MODEL_MP05: >>> - self.assertAlmostEqual(value, >>> 4909.88110195, 3) >>> - elif model == MODEL_NS_R1RHO_2SITE: >>> - self.assertAlmostEqual(value, >>> 5610.20221435, 2) >>> - >>> - elif param == 'chi2': >>> - if model == MODEL_NOREX: >>> - self.assertAlmostEqual(value, >>> 3363.95829122) >>> - elif model == MODEL_DPL94: >>> - self.assertAlmostEqual(value, >>> 710.24767560) >>> - elif model == MODEL_TP02: >>> - self.assertAlmostEqual(value, >>> 114.47142772) >>> - elif model == MODEL_TAP03: >>> - self.assertAlmostEqual(value, >>> 114.27987534) >>> - elif model == MODEL_MP05: >>> - self.assertAlmostEqual(value, >>> 114.28002272) >>> - elif model == MODEL_NS_R1RHO_2SITE: >>> - self.assertAlmostEqual(value, >>> 134.14368365) >>> + if do_assert: >>> + # Compare values. >>> + if spin_id == ':52@N': >>> + if param == 'phi_ex': >>> + if model == MODEL_DPL94: >>> + self.assertAlmostEqual(value, >>> 0.07561937) >>> + >>> + elif param == 'pA': >>> + if model == MODEL_TP02: >>> + self.assertAlmostEqual(value, >>> 0.88807487) >>> + elif model == MODEL_TAP03: >>> + self.assertAlmostEqual(value, >>> 0.88809318) >>> + elif model == MODEL_MP05: >>> + self.assertAlmostEqual(value, >>> 0.88809321) >>> + elif model == MODEL_NS_R1RHO_2SITE: >>> + self.assertAlmostEqual(value, >>> 0.94496541, 6) >>> + >>> + elif param == 'dw': >>> + if model == MODEL_TP02: >>> + self.assertAlmostEqual(value, >>> 1.08765638, 6) >>> + elif model == MODEL_TAP03: >>> + self.assertAlmostEqual(value, >>> 1.08726698, 6) >>> + elif model == MODEL_MP05: >>> + self.assertAlmostEqual(value, >>> 1.08726706, 6) >>> + elif model == MODEL_NS_R1RHO_2SITE: >>> + self.assertAlmostEqual(value, >>> 1.55833321, 5) >>> + >>> + elif param == 'kex': >>> + if model == MODEL_DPL94: >>> + self.assertAlmostEqual(value, >>> 4419.03917195, 2) >>> + elif model == MODEL_TP02: >>> + self.assertAlmostEqual(value, >>> 4904.70144883, 3) >>> + elif model == MODEL_TAP03: >>> + self.assertAlmostEqual(value, >>> 4909.86877150, 3) >>> + elif model == MODEL_MP05: >>> + self.assertAlmostEqual(value, >>> 4909.88110195, 3) >>> + elif model == MODEL_NS_R1RHO_2SITE: >>> + self.assertAlmostEqual(value, >>> 5610.20221435, 2) >>> + >>> + elif param == 'chi2': >>> + if model == MODEL_NOREX: >>> + self.assertAlmostEqual(value, >>> 3363.95829122) >>> + elif model == MODEL_DPL94: >>> + self.assertAlmostEqual(value, >>> 710.24767560) >>> + elif model == MODEL_TP02: >>> + self.assertAlmostEqual(value, >>> 114.47142772) >>> + elif model == MODEL_TAP03: >>> + self.assertAlmostEqual(value, >>> 114.27987534) >>> + elif model == MODEL_MP05: >>> + self.assertAlmostEqual(value, >>> 114.28002272) >>> + elif model == MODEL_NS_R1RHO_2SITE: >>> + self.assertAlmostEqual(value, >>> 134.14368365) >>> >>> >>> # Print the final pipe. >>> @@ -7524,7 +7529,7 @@ >>> >>> # Assign the split of the line. >>> mol_name, res_num, res_name, spin_num, spin_name, val, >>> sd_error = line_split >>> - print mol_name, res_num, res_name, spin_num, spin_name, >>> val, sd_error >>> + print(mol_name, res_num, res_name, spin_num, spin_name, >>> val, sd_error) >>> >>> if res_num == '52': >>> # Assert that the value is not None. >>> >>> >>> _______________________________________________ >>> 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
