r25781 - /trunk/pipe_control/error_analysis.py -- September 12, 2014 - 13:48

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):