This is a good test - bootstrapping results should be similar to the covariance matrix and MC simulations. However you should note that statisticians do not consider this technique to be a valid error estimate.
Regards, Edward On 31 August 2014 17:26, <[email protected]> wrote: > Author: tlinnet > Date: Sun Aug 31 17:26:48 2014 > New Revision: 25485 > > URL: http://svn.gna.org/viewcvs/relax?rev=25485&view=rev > Log: > Modified systemtest Relax_disp.verify_estimate_r2eff_err_compare_mc to > include boot strapping method. > > This shows there is excellent agreement between boot-strapping and estimation > of errors from Co-variance, while > relax Monte-Carlo simulations are very much different. > > Boot strapping is the "-2": > > -2 0.070 0.085 0.087 0.095 0.086 0.076 0.087 0.072 0.069 0.077 0.025 > 0.035 0.018 0.015 sum= 0.897 > -1 0.069 0.081 0.085 0.092 0.085 0.074 0.083 0.069 0.066 0.074 0.025 > 0.035 0.018 0.016 sum= 0.874 > 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 > 400 0.034 0.043 0.044 0.049 0.046 0.037 0.042 0.035 0.031 0.039 0.014 > 0.018 0.009 0.008 sum= 0.450 > 800 0.032 0.040 0.041 0.046 0.042 0.036 0.040 0.034 0.034 0.037 0.013 > 0.018 0.009 0.008 sum= 0.431 > 1200 0.033 0.041 0.042 0.046 0.043 0.037 0.042 0.036 0.034 0.038 0.012 > 0.018 0.009 0.008 sum= 0.439 > 1600 0.036 0.041 0.042 0.047 0.043 0.038 0.042 0.035 0.035 0.037 0.013 > 0.018 0.009 0.008 sum= 0.443 > 2000 0.034 0.042 0.042 0.046 0.042 0.036 0.043 0.035 0.034 0.037 0.013 > 0.017 0.009 0.008 sum= 0.438 > > 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=25485&r1=25484&r2=25485&view=diff > ============================================================================== > --- trunk/test_suite/system_tests/relax_disp.py (original) > +++ trunk/test_suite/system_tests/relax_disp.py Sun Aug 31 17:26:48 2014 > @@ -24,6 +24,8 @@ > # Python module imports. > from os import F_OK, access, getcwd, path, sep > from numpy import array, asarray, exp, median, inf, log, save, std, sum, > zeros > +from minfx.generic import generic_minimise > +from random import gauss > import re, math > from tempfile import mkdtemp, NamedTemporaryFile > > @@ -3125,9 +3127,7 @@ > # Do boot strapping ? > do_boot = True > if do_boot: > - from minfx.generic import generic_minimise > - from random import gauss > - min_algor = 'BFGS' > + min_algor = 'Newton' > min_options = () > sim_boot = 2000 > scaling_list = [1.0, 1.0] > @@ -3232,7 +3232,7 @@ > x0 = [r2eff, i0] > setup(num_params=len(x0), num_times=len(times), > values=I_err, sd=errors, relax_times=times, scaling_matrix=scaling_list) > > - params_minfx_sim_j, chi2_minfx_sim_j, iter_count, > f_count, g_count, h_count, warning = generic_minimise(func=func, dfunc=dfunc, > args=(), x0=x0, min_algor=min_algor, min_options=min_options, > full_output=True, print_flag=0) > + params_minfx_sim_j, chi2_minfx_sim_j, iter_count, > f_count, g_count, h_count, warning = generic_minimise(func=func, dfunc=dfunc, > d2func=d2func, args=(), x0=x0, min_algor=min_algor, min_options=min_options, > full_output=True, print_flag=0) > R_m_sim_j, I0_m_sim_j = params_minfx_sim_j > R_m_sim_l.append(R_m_sim_j) > I0_m_sim_l.append(I0_m_sim_j) > @@ -7995,6 +7995,9 @@ > # Set algorithm. > min_algor = 'Newton' > constraints = False > + min_options = () > + sim_boot = 2000 > + scaling_list = [1.0, 1.0] > > # Minimise. > self.interpreter.minimise.execute(min_algor=min_algor, > constraints=constraints, verbosity=1) > @@ -8008,10 +8011,7 @@ > if hasattr(cur_spin, err_attr): > delattr(cur_spin, err_attr) > > - # Estimate R2eff errors. > - self.interpreter.relax_disp.r2eff_err_estimate(chi2_jacobian=False) > - > - # Collect the estimation data. > + # Collect the estimation data from boot. > my_dic = {} > param_key_list = [] > est_keys = [] > @@ -8038,6 +8038,69 @@ > > # Append key. > param_key_list.append(param_key) > + > + # Add key to dic. > + my_dic[spin_id][est_key][param_key] = {} > + > + values = [] > + errors = [] > + times = [] > + for time in loop_time(exp_type=exp_type, frq=frq, > offset=offset, point=point): > + values.append(average_intensity(spin=cur_spin, > exp_type=exp_type, frq=frq, offset=offset, point=point, time=time)) > + errors.append(average_intensity(spin=cur_spin, > exp_type=exp_type, frq=frq, offset=offset, point=point, time=time, > error=True)) > + times.append(time) > + > + # Convert to numpy array. > + values = asarray(values) > + errors = asarray(errors) > + times = asarray(times) > + > + r2eff = getattr(cur_spin, 'r2eff')[param_key] > + i0 = getattr(cur_spin, 'i0')[param_key] > + > + R_m_sim_l = [] > + I0_m_sim_l = [] > + for j in range(sim_boot): > + if j in range(0, 100000, 100): > + print("Simulation %i"%j) > + # Start minimisation. > + > + # Produce errors > + I_err = [] > + for j, error in enumerate(errors): > + I_error = gauss(values[j], error) > + I_err.append(I_error) > + # Convert to numpy array. > + I_err = asarray(I_err) > + > + x0 = [r2eff, i0] > + setup(num_params=len(x0), num_times=len(times), > values=I_err, sd=errors, relax_times=times, scaling_matrix=scaling_list) > + > + params_minfx_sim_j, chi2_minfx_sim_j, iter_count, > f_count, g_count, h_count, warning = generic_minimise(func=func, dfunc=dfunc, > d2func=d2func, args=(), x0=x0, min_algor=min_algor, min_options=min_options, > full_output=True, print_flag=0) > + R_m_sim_j, I0_m_sim_j = params_minfx_sim_j > + R_m_sim_l.append(R_m_sim_j) > + I0_m_sim_l.append(I0_m_sim_j) > + > + # Get stats on distribution. > + sigma_R_sim = std(asarray(R_m_sim_l), ddof=1) > + sigma_I0_sim = std(asarray(I0_m_sim_l), ddof=1) > + my_dic[spin_id][est_key][param_key]['r2eff_err'] = > sigma_R_sim > + my_dic[spin_id][est_key][param_key]['i0_err'] = sigma_I0_sim > + > + # Estimate R2eff errors. > + self.interpreter.relax_disp.r2eff_err_estimate() > + > + est_key = '-1' > + est_keys.append(est_key) > + > + # Collect data. > + for cur_spin, mol_name, resi, resn, spin_id in > spin_loop(full_info=True, return_id=True, skip_desel=True): > + # Add key for estimate. > + my_dic[spin_id][est_key] = {} > + > + for exp_type, frq, offset, point, ei, mi, oi, di in > loop_exp_frq_offset_point(return_indices=True): > + # Generate the param_key. > + param_key = return_param_key_from_data(exp_type=exp_type, > frq=frq, offset=offset, point=point) > > # Add key to dic. > my_dic[spin_id][est_key][param_key] = {} > @@ -8054,38 +8117,8 @@ > my_dic[spin_id][est_key][param_key][err_attr] = > get_err_attr > > > - # Estimate R2eff errors from Chi2 Jacobian. > - self.interpreter.relax_disp.r2eff_err_estimate(chi2_jacobian=True) > - > - est_key = '-1' > - est_keys.append(est_key) > - > - # Collect data. > - for cur_spin, mol_name, resi, resn, spin_id in > spin_loop(full_info=True, return_id=True, skip_desel=True): > - # Add key for estimate. > - my_dic[spin_id][est_key] = {} > - > - for exp_type, frq, offset, point, ei, mi, oi, di in > loop_exp_frq_offset_point(return_indices=True): > - # Generate the param_key. > - param_key = return_param_key_from_data(exp_type=exp_type, > frq=frq, offset=offset, point=point) > - > - # Add key to dic. > - my_dic[spin_id][est_key][param_key] = {} > - > - # Get the value. > - # Loop over err attributes. > - for err_attr in err_attr_list: > - if hasattr(cur_spin, err_attr): > - get_err_attr = getattr(cur_spin, err_attr)[param_key] > - else: > - get_err_attr = 0.0 > - > - # Save to dic. > - my_dic[spin_id][est_key][param_key][err_attr] = > get_err_attr > - > - > # Make Carlo Simulations number > - mc_number_list = range(0, 50, 10) > + mc_number_list = range(0, 2400, 400) > > sim_attr_list = ['chi2_sim', 'f_count_sim', 'g_count_sim', > 'h_count_sim', 'i0_sim', 'iter_sim', 'peak_intensity_sim', 'r2eff_sim', > 'select_sim', 'warning_sim'] > > > > _______________________________________________ > 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
