Author: tlinnet
Date: Tue Sep 2 03:06:02 2014
New Revision: 25526
URL: http://svn.gna.org/viewcvs/relax?rev=25526&view=rev
Log:
Tried to implement bootstrap method.
But that goes horrible wrong.
The numbers are now:
param: r2a_err, with err: 0.48131663, compared to MC: 0.49354342, with 600
boot 0.00000000
param: dw_err, with err: 117.99734537, compared to MC: 0.32813703, with 600
boot 0.00000000
param: k_AB_err, with err: 0.64606994, compared to MC: 0.65384783, with 600
boot 0.00000000
So, something is right, and something is really wrong.
It is clear, that the reported error estimation scales with the error.
So, either one should start with simulations from peak intensity and then
forward.
Or, maybe the use of errors are wrong.
There are some rules for how errors can be added together.
If they are a product, then the fractional errors can be summed together.
But, I have to look into this.
task #7824(https://gna.org/task/index.php?7824): Model parameter ERROR
estimation from Jacobian and Co-variance matrix of dispersion models.
Modified:
branches/est_par_error/test_suite/system_tests/relax_disp.py
Modified: branches/est_par_error/test_suite/system_tests/relax_disp.py
URL:
http://svn.gna.org/viewcvs/relax/branches/est_par_error/test_suite/system_tests/relax_disp.py?rev=25526&r1=25525&r2=25526&view=diff
==============================================================================
--- branches/est_par_error/test_suite/system_tests/relax_disp.py
(original)
+++ branches/est_par_error/test_suite/system_tests/relax_disp.py Tue
Sep 2 03:06:02 2014
@@ -41,12 +41,13 @@
from pipe_control.minimise import assemble_scaling_matrix
from specific_analyses.relax_disp.checks import check_missing_r1
from specific_analyses.relax_disp.estimate_r2eff import estimate_par_err,
estimate_r2eff
-from specific_analyses.relax_disp.data import average_intensity,
check_intensity_errors, generate_r20_key, get_curve_type,
has_exponential_exp_type, has_r1rho_exp_type, loop_exp_frq,
loop_exp_frq_offset, loop_exp_frq_offset_point,
loop_exp_frq_offset_point_time, loop_time, return_grace_file_name_ini,
return_param_key_from_data, spin_ids_to_containers
+from specific_analyses.relax_disp.data import average_intensity,
check_intensity_errors, generate_r20_key, get_curve_type,
has_exponential_exp_type, has_r1rho_exp_type, is_r1_optimised, loop_exp_frq,
loop_exp_frq_offset, loop_exp_frq_offset_point,
loop_exp_frq_offset_point_time, loop_point, loop_time, return_cpmg_frqs,
return_grace_file_name_ini, return_param_key_from_data,
spin_ids_to_containers, return_offset_data, return_r1_data,
return_r1_err_data, return_r2eff_arrays, return_spin_lock_nu1
from specific_analyses.relax_disp.data import INTERPOLATE_DISP,
INTERPOLATE_OFFSET, X_AXIS_DISP, X_AXIS_W_EFF, X_AXIS_THETA, Y_AXIS_R2_R1RHO,
Y_AXIS_R2_EFF
from specific_analyses.relax_disp.model import models_info, nesting_param
-from specific_analyses.relax_disp.parameters import linear_constraints
+from specific_analyses.relax_disp.parameters import assemble_param_vector,
linear_constraints, param_num
from specific_analyses.relax_disp.variables import EXP_TYPE_CPMG_DQ,
EXP_TYPE_CPMG_MQ, EXP_TYPE_CPMG_PROTON_MQ, EXP_TYPE_CPMG_PROTON_SQ,
EXP_TYPE_CPMG_SQ, EXP_TYPE_CPMG_ZQ, EXP_TYPE_LIST, EXP_TYPE_R1RHO,
MODEL_B14_FULL, MODEL_CR72, MODEL_CR72_FULL, MODEL_DPL94, MODEL_IT99,
MODEL_LIST_ANALYTIC_CPMG, MODEL_LIST_FULL, MODEL_LIST_NUMERIC_CPMG,
MODEL_LM63, MODEL_M61, MODEL_M61B, MODEL_MP05, MODEL_NOREX,
MODEL_NS_CPMG_2SITE_3D_FULL, MODEL_NS_CPMG_2SITE_EXPANDED,
MODEL_NS_CPMG_2SITE_STAR_FULL, MODEL_NS_R1RHO_2SITE, MODEL_NS_R1RHO_3SITE,
MODEL_NS_R1RHO_3SITE_LINEAR, MODEL_PARAMS, MODEL_R2EFF, MODEL_TP02,
MODEL_TAP03, PARAMS_R20
from status import Status; status = Status()
+from target_functions.relax_disp import Dispersion
from test_suite.system_tests.base_classes import SystemTestCase
# C modules.
@@ -7377,6 +7378,38 @@
my_dic = {}
spin_id_list = []
+ # Set algorithm.
+ min_algor = 'simplex'
+ constraints = True
+ if constraints:
+ min_options = ('%s'%(min_algor),)
+ min_algor = 'Log barrier'
+ #min_algor = 'Method of Multipliers'
+ scaling_matrix = assemble_scaling_matrix(scaling=True)
+
+ # Collect spins
+ all_spin_ids = []
+ for cur_spin, mol_name, resi, resn, spin_id in
spin_loop(full_info=True, return_id=True, skip_desel=True):
+ all_spin_ids.append(spin_id)
+
+ spins = spin_ids_to_containers(all_spin_ids[:1])
+
+ # Get constraints
+ A, b = linear_constraints(spins=spins,
scaling_matrix=scaling_matrix[0])
+ else:
+ min_options = ()
+ A, b = None, None
+
+
+ sim_boot = 2
+
+ # Number of spectrometer fields.
+ fields = [None]
+ field_count = 1
+ if hasattr(cdp, 'spectrometer_frq_count'):
+ fields = cdp.spectrometer_frq_list
+ field_count = cdp.spectrometer_frq_count
+
# Get the data.
for cur_spin, mol_name, resi, resn, spin_id in
spin_loop(full_info=True, return_id=True, skip_desel=True):
# Add key to dic.
@@ -7392,6 +7425,10 @@
# Generate spin string.
spin_string = generate_spin_string(spin=cur_spin,
mol_name=mol_name, res_num=resi, res_name=resn)
+ # Get the param_vector
+ param_vector = assemble_param_vector(spins=[cur_spin])
+ model = cur_spin.model
+
param_key_list = []
for exp_type, frq, offset, ei, mi, oi, in
loop_exp_frq_offset(return_indices=True):
# Get the parameter key.
@@ -7418,6 +7455,73 @@
param_err_val = deepcopy(getattr(cur_spin,
param_err))
my_dic[spin_id][param_key][param_err] = param_err_val
+ # Start boot strap.
+ values, errors, missing, frqs, frqs_H, exp_types,
relax_times = return_r2eff_arrays(spins=[cur_spin], spin_ids=[spin_id],
fields=fields, field_count=field_count)
+ # The offset data.
+ offsets, spin_lock_fields_inter, chemical_shifts,
tilt_angles, Delta_omega, w_eff = return_offset_data(spins=[cur_spin],
spin_ids=[spin_id], field_count=field_count)
+
+ # r1 data.
+ r1 = return_r1_data(spins=[cur_spin], spin_ids=[spin_id],
field_count=field_count)
+ r1_err = return_r1_err_data(spins=[cur_spin],
spin_ids=[spin_id], field_count=field_count)
+ r1_fit = is_r1_optimised(model)
+ model_param_num = param_num(spins=[cur_spin])
+
+ dispersion_points = cdp.dispersion_points
+ cpmg_frqs = return_cpmg_frqs(ref_flag=False)
+ spin_lock_nu1 = return_spin_lock_nu1(ref_flag=False)
+ num_spins=1
+
+ values_err = deepcopy(values)
+ si = 0
+ r2a_sim_l = []
+ dw_sim_l = []
+ k_AB_sim_l = []
+
+ for j in range(sim_boot):
+ if j in range(0, 1000, 10):
+ print("Simulation %i"%j)
+ # Start minimisation.
+
+ # Produce errors
+ r2_err = []
+ for j, error in enumerate(errors[ei][si][mi][oi]):
+ r2_error = gauss(values[ei][si][mi][oi][j], error)
+ #print values[ei][si][mi][oi][j], r2_error, error
+ values_err[ei][si][mi][oi][j] = deepcopy(r2_error)
+
+ # Init the Dispersion clas with data, so we can call
functions in it.
+ tfunc = Dispersion(model=model,
num_params=model_param_num, num_spins=num_spins, num_frq=field_count,
exp_types=exp_types, values=values_err, errors=errors, missing=missing,
frqs=frqs, frqs_H=frqs_H, cpmg_frqs=cpmg_frqs, spin_lock_nu1=spin_lock_nu1,
chemical_shifts=chemical_shifts, offset=offsets, tilt_angles=tilt_angles,
r1=r1, relax_times=relax_times, r1_fit=r1_fit)
+ results = generic_minimise(func=tfunc.func, args=(),
x0=param_vector, min_algor=min_algor, min_options=min_options, A=A, b=b,
full_output=True, print_flag=0)
+ param_vector, chi2, iter_count, f_count, g_count,
h_count, warning = results
+ print chi2, warning
+
+ # Loop over params.
+ for i, param in enumerate(params):
+ # If param in PARAMS_R20, values are stored in with
parameter key.
+ if param in PARAMS_R20:
+ # Get the error.
+ param_val = deepcopy(getattr(cur_spin,
param)[param_key])
+ r2a_sim_l.append(param_val)
+
+ else:
+ # Get the error.
+ param_val = deepcopy(getattr(cur_spin, param))
+ if param == 'dw':
+ dw_sim_l.append(param_val)
+ elif param == 'k_AB':
+ k_AB_sim_l.append(param_val)
+
+ # Get stats on distribution.
+ sigma_r2a_sim = std(asarray(r2a_sim_l), ddof=1)
+ sigma_dw_sim = std(asarray(dw_sim_l), ddof=1)
+ sigma_k_AB_sim = std(asarray(k_AB_sim_l), ddof=1)
+ my_dic[spin_id][param_key]['sigma_r2a_sim'] = sigma_r2a_sim
+ my_dic[spin_id][param_key]['sigma_dw_sim'] = sigma_dw_sim
+ my_dic[spin_id][param_key]['sigma_k_AB_sim'] = sigma_k_AB_sim
+ my_dic[spin_id][param_key]['sigma_r2a_arr'] =
asarray(r2a_sim_l)
+ my_dic[spin_id][param_key]['sigma_dw_arr'] =
asarray(dw_sim_l)
+ my_dic[spin_id][param_key]['sigma_k_AB_arr'] =
asarray(k_AB_sim_l)
+
my_dic[spin_id]['param_key_list'] = param_key_list
# After Monte-Carlo.
@@ -7456,11 +7560,14 @@
param_key_list = my_dic[spin_id]['param_key_list']
for param_key in param_key_list:
params_err = my_dic[spin_id]['params_err']
+ params = my_dic[spin_id]['params']
# Loop over params.
for i, param_err in enumerate(params_err):
+ param = params[i]
param_err_val = my_dic[spin_id][param_key][param_err]
param_err_val_mc =
my_dic[spin_id][param_key][param_err+'MC']
- print("param: %s, with err: %3.8f, compared to MC:
%3.8f" % (param_err, param_err_val, param_err_val_mc) )
+ param_err_val_sim =
my_dic[spin_id][param_key]['sigma_%s_sim'%param]
+ print("param: %s, with err: %3.8f, compared to MC:
%3.8f, with %i boot %3.8f" % (param_err, param_err_val, param_err_val_mc,
sim_boot, param_err_val_sim) )
def test_tp02_data_to_ns_r1rho_2site(self, model=None):
r25526 - /branches/est_par_error/test_suite/system_tests/relax_disp.py