Hi Troels, I'm not sure if I can answer your question about the rows/columns of the Jacobian in the code here.
Regards, Edward On 1 September 2014 23:25, <[email protected]> wrote: > Author: tlinnet > Date: Mon Sep 1 23:25:36 2014 > New Revision: 25519 > > URL: http://svn.gna.org/viewcvs/relax?rev=25519&view=rev > Log: > Further extended the script for analysing errors with Jacobian. > > Now loops over the dimension. > > 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/specific_analyses/relax_disp/estimate_r2eff.py > > Modified: > branches/est_par_error/specific_analyses/relax_disp/estimate_r2eff.py > URL: > http://svn.gna.org/viewcvs/relax/branches/est_par_error/specific_analyses/relax_disp/estimate_r2eff.py?rev=25519&r1=25518&r2=25519&view=diff > ============================================================================== > --- branches/est_par_error/specific_analyses/relax_disp/estimate_r2eff.py > (original) > +++ branches/est_par_error/specific_analyses/relax_disp/estimate_r2eff.py > Mon Sep 1 23:25:36 2014 > @@ -38,9 +38,9 @@ > from pipe_control.minimise import assemble_scaling_matrix > from pipe_control.mol_res_spin import generate_spin_string, spin_loop > from specific_analyses.relax_disp.checks import check_model_type > -from specific_analyses.relax_disp.data import average_intensity, > is_r1_optimised, loop_exp_frq_offset_point, loop_time, return_cpmg_frqs, > return_offset_data, return_param_key_from_data, return_r1_data, > return_r1_err_data, return_r2eff_arrays, return_spin_lock_nu1 > +from specific_analyses.relax_disp.data import average_intensity, > generate_r20_key, is_r1_optimised, loop_exp_frq_offset, > loop_exp_frq_offset_point, loop_time, return_cpmg_frqs, return_offset_data, > return_param_key_from_data, return_r1_data, return_r1_err_data, > return_r2eff_arrays, return_spin_lock_nu1 > from specific_analyses.relax_disp.parameters import assemble_param_vector, > disassemble_param_vector, param_num > -from specific_analyses.relax_disp.variables import MODEL_CR72, MODEL_R2EFF, > MODEL_TSMFK01 > +from specific_analyses.relax_disp.variables import MODEL_CR72, MODEL_R2EFF, > MODEL_TSMFK01, PARAMS_R20 > from target_functions.chi2 import chi2_rankN, dchi2 > from target_functions.relax_disp import Dispersion > > @@ -243,8 +243,53 @@ > jacobian = tfunc.jacobian(param_vector) > weights = 1. / tfunc.errors**2 > > - # Get the co-variance > - pcov = multifit_covar(J=jacobian, weights=weights) > + # Get the shape of the data. > + NJ, NE, NS, NM, NO, ND = jacobian.shape > + if NS != 1: > + raise RelaxError("The number of spins does not fit.") > + > + # Get the parameters fitted in the model. > + params = cur_spin.params > + > + # Set the spin index to 0. > + si = 0 > + # Loop over the data. > + for exp_type, frq, offset, ei, mi, oi in > loop_exp_frq_offset(return_indices=True): > + param_key = generate_r20_key(exp_type=exp_type, frq=frq) > + > + # Extract weights. > + cur_weights = weights[ei, si, mi, oi] > + > + # Extract every column/row from the first to last columns. Is > this correct? > + cur_jacobian = jacobian[0:NJ:1, ei, si, mi, oi] > + > + # Get the co-variance > + pcov = multifit_covar(J=cur_jacobian, weights=cur_weights) > + > + # To compute one standard deviation errors on the parameters, > take the square root of the diagonal covariance. > + param_vector_error = sqrt(diag(pcov)) > + > + # Loop over params. > + for i, param in enumerate(params): > + # Set the param error name > + param_err = param + '_err' > + > + # If param in PARAMS_R20, values are stored in with > parameter key. > + if param in PARAMS_R20: > + # Copy parameter attribute to error attribute, if not in > spin Class. > + if not hasattr(cur_spin, param_err): > + setattr(cur_spin, param_err, > deepcopy(getattr(cur_spin, param))) > + > + # Set error. > + getattr(cur_spin, param_err)[param_key] = > deepcopy(param_vector_error[i]) > + > + else: > + # Copy parameter attribute to error attribute, if not in > spin Class. > + if not hasattr(cur_spin, param_err): > + setattr(cur_spin, param_err, > deepcopy(getattr(cur_spin, param))) > + > + # Set error. > + setattr(cur_spin, param_err, param_vector_error[i]) > > > #### This class is only for testing. > > > _______________________________________________ > 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
