Changeset 4efc00a in tol

Timestamp:

May 11, 2010, 2:58:02 PM (16 years ago)

Author:

Víctor de Buen <vdebuen@…>

Branches:

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

Children:

a868db9

Parents:

5633d43

Message:

multivariate linear model with unknown covariance matrix and constraining linear equalities
refs #902

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

Location:

tol_pkg/BysVecLinReg

Files:

• Multivariate.tol (modified) (3 diffs)
• test/gen_multivariant_linear_regression.tol (modified) (2 diffs)
• test/test_0003/test.tol (added)
• tools.tol (modified) (1 diff)

tol_pkg/BysVecLinReg/Multivariate.tol

@@ -73 +73 @@
   //Choleski decomposition of VxIm as sparse explicit((d*m)x(d*m))
   VMatrix _.VxIm.Le = Rand(0,0,0,0);
- 
+
+  ////////////////////////////////////////////////////////////////////////////
+  // Optional linear constraining equalities
+  // 
+  //   C * vec(A) = c
+  // 
+  //Number of linear constraining equalities
+  Real _.r.C = 0;
+  //Coefficient matrix of linear constraining equalities 
+  VMatrix _.C = Rand(0,0,0,0);
+  //Right vector of linear constraining equalities 
+  VMatrix _.c = Rand(0,0,0,0);
+  ////////////////////////////////////////////////////////////////////////////
+
   ////////////////////////////////////////////////////////////////////////////
   Real check.dim(Text method, VMatrix M, Real nRow, Real nCol)
@@ -155 +168 @@
 
   ////////////////////////////////////////////////////////////////////////////
+  Real set.linearConstrainingEquations(VMatrix C, VMatrix c)
+  ////////////////////////////////////////////////////////////////////////////
+  {
+    Real r = VRows(c);
+    If(Or(!check.dim("set.linearConstrainingEquations",C,r,_.m*_.d),
+          !check.dim("set.linearConstrainingEquations",c,r,1      )),False,
+    {
+      Real _.r.C := r;
+      VMatrix _.C := C; 
+      VMatrix _.c := c; 
+      True
+    })
+  };
+
+  ////////////////////////////////////////////////////////////////////////////
   Real change.X(VMatrix X)
   ////////////////////////////////////////////////////////////////////////////
@@ -251 +279 @@
 
   ////////////////////////////////////////////////////////////////////////////
+  VMatrix gibbs.draw.vec.A.std(Real numSim)
+  ////////////////////////////////////////////////////////////////////////////
+  {
+    If(!_.r.C, Gaussian(_.d*_.m,1,0,1),
+    {
+      VMatrix C = Tra(_.IdxXLi*Tra(_.C))*_.VxIm.Le;
+      VMatrix c = _.c - _.C * _.vec.nu;
+      eqConstrainedStdGaussian(C,c)
+    })
+  };
+
+  ////////////////////////////////////////////////////////////////////////////
   VMatrix gibbs.draw.vec.A(Real numSim)
   ////////////////////////////////////////////////////////////////////////////
   {
-    VMatrix e = Gaussian(_.d*_.m,1,0,1);
+    VMatrix e = gibbs.draw.vec.A.std(numSim);
     VMatrix e.1 = _.VxIm.Le*e;
     VMatrix e.2 = Tra(Tra(e.1)*_.IdxXLi);

tol_pkg/BysVecLinReg/test/gen_multivariant_linear_regression.tol

@@ -10 +10 @@
 WriteLn("Current Random Seed = "<<rndSeed);
 
-Real d = IntRand(2,10);
-Real m = IntRand(2,10);
-Real N = IntRand(m*10,m*100);
+Real If(!ObjectExist("Real","r.eq"), r.eq = 0, r.eq);
+
+Real d = IntRand(r.eq+2,r.eq+10);
+Real m = IntRand(r.eq+2,r.eq+10);
+Real N = IntRand(m*5,m*30);
 Real sampleLength = 5000;
 
@@ -23 +25 @@
 VMatrix Y = A*X+E; 
 
+VMatrix C = If(!r.eq,Rand(0,d*m,0,0), 
+  Round(Rand(r.eq,d*m,-2.5,2.5,"Cholmod.R.Sparse",r.eq*2)));
+VMatrix c = If(!r.eq,Rand(0,1,0,0),
+  C*BysVecLinReg::vmatrix2vec(A));

tol_pkg/BysVecLinReg/tools.tol

@@ -43 +43 @@
   VMatrix c.2 = c-c.1;
   VMatrix e.2 = MinimumResidualsSolve(C.2,c.2);
-//VMatrix C.2*e.2-c.2
-  e.2
+  VMatrix err = C.2*e.2-c.2;
+  Real errMax = VMatMax(Abs(err));
+  If((n^3)+errMax^2 != n^3,
+    WriteLn("[BysVecLinReg::eqConstrainedStdGaussian] errMax="<<errMax,"W"));
+  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))
 };