Changeset e9469de4 in tol

Timestamp:

Mar 27, 2009, 10:32:20 AM (18 years ago)

Author:

Víctor de Buen <vdebuen@…>

Branches:

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

Children:

44ae460

Parents:

76ecf8e

Message:

Preparing for Metropolis-Hastings method explained in #664

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

Location:

tol/stdlib/math/stat/models/bayesian/bysMcmc/bsr/gibbs

Files:

• _basicMaster.tol (modified) (1 diff)
• _inputMissingBlock.tol (modified) (1 diff)
• _outputMissingBlock.tol (modified) (1 diff)
• _stdLinearBlock.tol (modified) (11 diffs)

tol/stdlib/math/stat/models/bayesian/bysMcmc/bsr/gibbs/_basicMaster.tol

@@ -264 +264 @@
   //WriteLn("TRACE"+_MID+" draw.sig.lin 8");
     VMatrix _.beta := If(Not(_.lin.blk::_.blockDef::_.numCol),Rand(0,1,0,0),
-      _.lin.blk::draw(numSim, Y01, X01));
+    {
+      Real _.lin.blk::solve(numSim, Y01, X01);
+      _.lin.blk::draw(numSim, Y01, X01)
+    });
   //WriteLn("TRACE"+_MID+" draw.sig.lin 9");
     _.beta << sig << arm << im << om

tol/stdlib/math/stat/models/bayesian/bysMcmc/bsr/gibbs/_inputMissingBlock.tol

@@ -210 +210 @@
     {
       NameBlock ls = buildLinearSystem(arima,si,beta,Z);
+      Real _.lin.blk::solve(numSim, ls::Y01, ls::X01);
       VMatrix _.blockDef::_.store := _.lin.blk::draw(numSim, ls::Y01, ls::X01);
       _.blockDef::_.store

tol/stdlib/math/stat/models/bayesian/bysMcmc/bsr/gibbs/_outputMissingBlock.tol

@@ -191 +191 @@
     {
       NameBlock ls = buildLinearSystem(arima,si,Z);
+      Real _.lin.blk::solve(numSim, ls::Y01, ls::X01);
       VMatrix _.blockDef::_.store := _.lin.blk::draw(numSim, ls::Y01, ls::X01);
       _.blockDef::_.store

tol/stdlib/math/stat/models/bayesian/bysMcmc/bsr/gibbs/_stdLinearBlock.tol

@@ -56 +56 @@
   //Matrix of linear system in last iteration
   VMatrix _old.X  = Rand(0,0,0,0); 
+  //
+  Set _.solution = Copy(Empty);
+
   Real _changed.X = True;
 
@@ -64 +67 @@
   //Current maximum of rounding error in solution
   Real _.cholSolMaxRoundErr     = 0;
+
+  //Stores the associated drawing method
+  Code _drawer = UnknownCode;
 
   ////////////////////////////////////////////////////////////////////////////
@@ -189 +195 @@
   ////////////////////////////////////////////////////////////////////////////
   //Solves the linear system corresponding to the standarized linear regression
-  NameBlock solve(Real numSim, VMatrix Y, VMatrix X)
+  Real solve(Real numSim, VMatrix Y, VMatrix X)
   ////////////////////////////////////////////////////////////////////////////
   {
@@ -221 +227 @@
     });
   //WriteLn("TRACE "+_MID+" solve 8");
-    NameBlock [[ XtY, nu, err ]]
+    Set _.solution := [[ numSim, Y, X, XtY, nu, err ]];
+    Real SetIndexByName(_.solution);
+    1
   };
 
@@ -227 +235 @@
   Text _.autodoc.member.draw.unconstrained =
   "Method of generation of a linear block in unconstrained case";
-  VMatrix draw.unconstrained(Real numSim, VMatrix Y, VMatrix X, NameBlock sr)
+  VMatrix draw.unconstrained(Real numSim)
   ////////////////////////////////////////////////////////////////////////////
   {
@@ -236 +244 @@
   //WriteLn("TRACE "+_MID+" draw.unconstrained 4 _.L="<<Matrix(VMat2Mat(_.L)));
   //WriteLn("TRACE "+_MID+" draw.unconstrained 4 b2="<<Matrix(VMat2Mat(Tra(b2))));
-    VMatrix  (_.blockDef::_.store) := sr::nu+b2;
+    VMatrix  (_.blockDef::_.store) := _.solution::nu+b2;
   //WriteLn("TRACE "+_MID+" draw.unconstrained 5");
     Copy(_.blockDef::_.store)
@@ -244 +252 @@
   Text _.autodoc.member.draw.constrained =
   "Method of generation of a linear block in constrained case";
-  VMatrix draw.constrained(Real numSim, VMatrix Y, VMatrix X, NameBlock sr)
-  ////////////////////////////////////////////////////////////////////////////
-  {
+  VMatrix draw.constrained(Real numSim)
+  ////////////////////////////////////////////////////////////////////////////
+  {
+    VMatrix nu = _.solution::nu;
   //WriteLn("TRACE "+_MID+" draw.constrained 1.1 Max(A*b0-a)="<<VMatMax(_.A*_.blockDef::_.store - a));
-    VMatrix d   = _.a - _.A*sr::nu;
+    VMatrix d   = _.a - _.A*nu;
     VMatrix If(_changed.X, _.D  := Tra(CholeskiSolve(_.L,_.At,"PtL")));
-    VMatrix z0  = Tra(Tra(_.blockDef::_.store-sr::nu)*_.Ls);
+    VMatrix z0  = Tra(Tra(_.blockDef::_.store-nu)*_.Ls);
   //WriteLn("TRACE "+_MID+" draw.constrained 1.2 Max(D*z0-d)="<<VMatMax(_.D*z0 - d));
     VMatrix z = TruncStdGaussian(_.D,d,z0,1,
@@ -265 +274 @@
       If(!VRows(_.beta0),0,
       {
-        VMatrix z0 := Tra(Tra(_.beta0-sr::nu)*_.Ls);
+        VMatrix z0 := Tra(Tra(_.beta0-nu)*_.Ls);
         VMatrix z  := TruncStdGaussian(_.D,d,z0,1,
                      _.blockDef::_.config::bsr.truncMNormal.gibbsNumIter);
@@ -272 +281 @@
       1
     });
-    VMatrix  (_.blockDef::_.store) := sr::nu + CholeskiSolve(_.L,z,"LtP");
+    VMatrix  (_.blockDef::_.store) := nu + CholeskiSolve(_.L,z,"LtP");
     Copy(_.blockDef::_.store)
   };
-
-  Code _drawer = draw.constrained;
 
   ////////////////////////////////////////////////////////////////////////////
@@ -286 +293 @@
     If(_.blockDef::enabled,
     {
-      NameBlock sr  = solve(numSim,Y,X);
-      VMatrix aux = _drawer(numSim, Y, X, sr);
+    //Real solve(numSim,Y,X);
+      VMatrix aux = _drawer(numSim);
     //WriteLn("TRACE "+_MID+" draw");
       Real _.blockDef::checkNonFiniteValues(numSim);
@@ -303 +310 @@
   ////////////////////////////////////////////////////////////////////////////
   {
-    NameBlock sr = solve(numSim,Y,X);
-  //WriteLn("TRACE "+_MID+" maxLikelihood sr::nu="<<Matrix(Tra(VMat2Mat(sr::nu))));
+    Real solve(numSim,Y,X);
     VMatrix aux = If(!_.R, 
     {
-      sr::nu
+      _.solution::nu
     },
     {
-      VMatrix d = _.a - _.A*sr::nu;
-      If(VMatMax(d)>=0, sr::nu, 
+      VMatrix d = _.a - _.A*_.solution::nu;
+      If(VMatMax(d)>=0, _.solution::nu, 
       {
         WriteLn("Cannot do maxLikelihood for constrained linear block "+
           _.blockDef::_.name,"W");
-        draw.constrained(numSim, Y, X, sr)
+        draw.constrained(numSim)
       })
     });