r25453 - in /trunk/test_suite/shared_data/curve_fitting/numeric_topology: covariance_matrix.log covariance_matrix.py -- August 29, 2014 - 20:18

Author: bugman
Date: Fri Aug 29 20:18:55 2014
New Revision: 25453

URL: http://svn.gna.org/viewcvs/relax?rev=25453&view=rev
Log:
Added a script and log for calculating the numerical covariance matrix for an 
exponential curve.

This uses the data at 
http://thread.gmane.org/gmane.science.nmr.relax.devel/6807/focus=6840 and
calculates the covariance matrix via the Jacobian calculated using the 
numdifftools.Jacobian object
construct and obtain the matrix, both at the minimum and at a point away from 
the minimum.  The
covariance is calculated as inv(J^T.W.J).


Added:
    
trunk/test_suite/shared_data/curve_fitting/numeric_topology/covariance_matrix.log
    
trunk/test_suite/shared_data/curve_fitting/numeric_topology/covariance_matrix.py
      - copied, changed from r25451, 
trunk/test_suite/shared_data/curve_fitting/numeric_topology/jacobian.py

Added: 
trunk/test_suite/shared_data/curve_fitting/numeric_topology/covariance_matrix.log
URL: 
http://svn.gna.org/viewcvs/relax/trunk/test_suite/shared_data/curve_fitting/numeric_topology/covariance_matrix.log?rev=25453&view=auto
==============================================================================
--- 
trunk/test_suite/shared_data/curve_fitting/numeric_topology/covariance_matrix.log
   (added)
+++ 
trunk/test_suite/shared_data/curve_fitting/numeric_topology/covariance_matrix.log
   Fri Aug 29 20:18:55 2014
@@ -0,0 +1,20 @@
+
+
+On-minimum:
+
+The covariance matrix at [1.0, 1000.0] is:
+[[  4.80710732e-04   7.51305845e-02]
+ [  7.51305845e-02   9.82126040e+01]]
+The parameter errors are therefore:
+    sigma_R:      0.02192511647039939796
+    sigma_I0:     9.91022724387324061013
+
+
+Off-minimum:
+
+The covariance matrix at [2.0, 500.0] is:
+[[  2.06611607e-02   1.92741253e-01]
+ [  1.92741253e-01   9.99664568e+01]]
+The parameter errors, which are rubbish as this is off-minimum, are 
therefore:
+    sigma_R:      0.14373990627293306566
+    sigma_I0:     9.99832270092464803213

Copied: 
trunk/test_suite/shared_data/curve_fitting/numeric_topology/covariance_matrix.py
 (from r25451, 
trunk/test_suite/shared_data/curve_fitting/numeric_topology/jacobian.py)
URL: 
http://svn.gna.org/viewcvs/relax/trunk/test_suite/shared_data/curve_fitting/numeric_topology/covariance_matrix.py?p2=trunk/test_suite/shared_data/curve_fitting/numeric_topology/covariance_matrix.py&p1=trunk/test_suite/shared_data/curve_fitting/numeric_topology/jacobian.py&r1=25451&r2=25453&rev=25453&view=diff
==============================================================================
--- trunk/test_suite/shared_data/curve_fitting/numeric_topology/jacobian.py   
  (original)
+++ 
trunk/test_suite/shared_data/curve_fitting/numeric_topology/covariance_matrix.py
    Fri Aug 29 20:18:55 2014
@@ -1,8 +1,10 @@
-# Script for numerically calculating the exponential curve gradient.
+# Script for numerically calculating the exponential curve covariance matrix 
via the Jacobian.
+# The equation used is:  covar = inv(J^T.W.J)
 
 # Python module imports.
-from math import exp
-from numpy import array, float64
+from math import exp, sqrt
+from numpy import array, dot, eye, float64, transpose
+from numpy.linalg import inv
 from numdifftools import Jacobian
 
 
@@ -36,17 +38,28 @@
 # The intensity errors.
 errors = [10.0, 10.0, 10.0, 10.0, 10.0]
 
+# The variance weighting matrix.
+W = eye(len(errors))/(10.0**2)
+
 # Set up the Jacobian function.
 jacobian = Jacobian(func)
 
-# The numeric Jacobian at the minimum.
+# The covariance matrix at the minimum.
 print("\n\nOn-minimum:\n")
-matrix = jacobian([R, I0])
-print("The Jacobian at %s is:\n%s" % ([R, I0], matrix))
+J = jacobian([R, I0])
+covar = inv(dot(transpose(J), dot(W, J)))
+print("The covariance matrix at %s is:\n%s" % ([R, I0], covar))
+print("The parameter errors are therefore:")
+print("    sigma_R:   %25.20f" % sqrt(covar[0, 0]))
+print("    sigma_I0:  %25.20f" % sqrt(covar[1, 1]))
 
-# The numeric Jacobian off the minimum.
+# The covariance matrix off the minimum.
 print("\n\nOff-minimum:\n")
 R = 2.0
 I0 = 500.0
-matrix = jacobian([R, I0])
-print("The Jacobian at %s is:\n%s" % ([R, I0], matrix))
+J = jacobian([R, I0])
+covar = inv(dot(transpose(J), dot(W, J)))
+print("The covariance matrix at %s is:\n%s" % ([R, I0], covar))
+print("The parameter errors, which are rubbish as this is off-minimum, are 
therefore:")
+print("    sigma_R:   %25.20f" % sqrt(covar[0, 0]))
+print("    sigma_I0:  %25.20f" % sqrt(covar[1, 1]))