Author: bugman Date: Fri Sep 12 13:48:58 2014 New Revision: 25781 URL: http://svn.gna.org/viewcvs/relax?rev=25781&view=rev Log: Implemented the pipe_control.error_analysis.covariance_matrix() function.
This follows from http://thread.gmane.org/gmane.science.nmr.relax.scm/23526/focus=7096. It will be used by a new error_analysis.covariance_matrix user function. And it calls the specific API methods model_loop(), covariance_matrix(), and set_error() and the relax library lib.statistics.multifit_covar() function do to most of the work. Modified: trunk/pipe_control/error_analysis.py Modified: trunk/pipe_control/error_analysis.py URL: http://svn.gna.org/viewcvs/relax/trunk/pipe_control/error_analysis.py?rev=25781&r1=25780&r2=25781&view=diff ============================================================================== --- trunk/pipe_control/error_analysis.py (original) +++ trunk/pipe_control/error_analysis.py Fri Sep 12 13:48:58 2014 @@ -31,6 +31,39 @@ from lib.errors import RelaxError from pipe_control import pipes from specific_analyses.api import return_api + + +def covariance_matrix(epsrel=0.0, verbosity=2): + """Estimate model parameter errors via the covariance matrix technique. + + Note that the covariance matrix error estimate is always of lower quality than Monte Carlo simulations. + + + @param epsrel: Any columns of R which satisfy |R_{kk}| <= epsrel |R_{11}| are considered linearly-dependent and are excluded from the covariance matrix, where the corresponding rows and columns of the covariance matrix are set to zero. + @type epsrel: float + @keyword verbosity: The amount of information to print. The higher the value, the greater the verbosity. + @type verbosity: int + """ + + # Test if the current data pipe exists. + pipes.test() + + # The specific analysis API object. + api = return_api() + + # Loop over the models. + for model_info in api.model_loop(): + # Get the Jacobian and weighting matrix. + jacobian, weights = api.covariance_matrix(verbosity=verbosity) + + # Calculate the covariance matrix. + pcov = statistics.multifit_covar(J=jacobian, weights=weights) + + # To compute one standard deviation errors on the parameters, take the square root of the diagonal covariance. + sd = sqrt(diag(pcov)) + + # Set the parameter error. + api.set_error(0, sd, model_info=model_info) def monte_carlo_create_data(method=None): _______________________________________________ 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
