Changeset 0c0e3bc in tol

Timestamp:

May 11, 2010, 4:31:09 PM (16 years ago)

Author:

Víctor de Buen <vdebuen@…>

Branches:

1941-docker-image-to-test, devel, master

Children:

ad7e8ca

Parents:

8d33ed4f

Message:

explaining eqConstrainedStdGaussian
refs #902

git-svn-id: ​https://www.tol-project.org/svn/tolp/trunk@2154 be4846ae-f051-483a-831f-b64e7d08527f

File:

• tol_pkg/BysVecLinReg/tools.tol (modified) (3 diffs)

tol_pkg/BysVecLinReg/tools.tol

@@ -47 +47 @@
 // 
 // 1. Get pivoting order p of columns of diagonal of C*C'
-// 3. Select p.1 as gaussian components of x as (n-r) last pivoted columns
-// 2. Select p.2 as dependent components of x as first r pivoted columns
+// 2. Select p.1 as gaussian components of x as (n-r) last pivoted columns
+// 3. Select p.2 as dependent components of x as first r pivoted columns
 // 4. Draw e.1 as standard independent normal 
 // 5. Rewrite constraining linear equalities conditioned to e.1
@@ -55 +55 @@
 //      e' = (e.1' e.2')
 //      C.1 * e.1 + C.2 * e.2 = c
-//      c.2 = c - C.1
+//      c.2 = c - C.1*e.1
 //      C.2 * e.2 = c.2
 // 
 // 6. Solve linear system to get e.2
-// 7. Merge e.1 and e.2 following the selected pivoting order to get the
+// 7. Checks if solution has the expected rounding error
+// 8. Merge e.1 and e.2 following the selected pivoting order to get the
 //    wanted draw of e
 // 
@@ -67 +68 @@
   Real r = VRows(C);
   Real n = VColumns(C);
-  VMatrix C2 = Constant(1,r,1) * (C$*C);
+  //Calculates diagonal of C*C' in a faster way
+  VMatrix C2 = Constant(1,r,1) * (C $* C);
+  //Builds pivoting order of columns of diagonal of C*C'
   Matrix p = ColumnPivots(VMat2Mat(C2));
+  //Select p.1 as gaussian components of x as (n-r) last pivoted columns
   Set p.1 = MatSet(Sub(p,1,r+1,1,n-r))[1];
+  //Select p.2 as dependent components of x as first r pivoted columns
   Set p.2 = MatSet(Sub(p,1,1,  1,r  ))[1];
-
+  //Draw e.1 as standard independent normal 
   VMatrix e.1 = Gaussian(n-r,1,0,1);
+  //Rewrite constraining linear equalities conditioned to e.1
   VMatrix C.1 = SubCol(C,p.1);
   VMatrix C.2 = SubCol(C,p.2);
-  VMatrix c.1 = C.1 * e.1;
-  VMatrix c.2 = c-c.1;
+  VMatrix c.2 = c-C.1 * e.1;
+  //Solve linear system to get e.2
   VMatrix e.2 = MinimumResidualsSolve(C.2,c.2);
+  //Checking method success
   VMatrix err = C.2*e.2-c.2;
   Real errMax = VMatMax(Abs(err));
-  If((n^3)+errMax^2 != n^3,
+  If((n*r^2)+errMax^2 != n*r^2,
     WriteLn("[BysVecLinReg::eqConstrainedStdGaussian] errMax="<<errMax,"W"));
+  //Merging gaussian and dependent components following the selected 
+  //pivoting order
   Triplet(SetCol(p.1<<p.2) | Constant(n,1,1) | VMat2Mat(e.1 << e.2), n, 1)
-//SetOfMatrix(Tra(p), Constant(n,1,1),VMat2Mat(e.1 << e.2))
 };