Changeset 40a7081 in tol

Timestamp:

Mar 27, 2009, 11:59:53 AM (18 years ago)

Author:

Víctor de Buen <vdebuen@…>

Branches:

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

Children:

96696c3

Parents:

44ae460

Message:

Preparing for Metropolis-Hastings method explained in #664

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

Location:

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

Files:

• _basicMaster.tol (modified) (8 diffs)
• _gibbs.tol (modified) (1 diff)
• _inputMissingBlock.tol (modified) (1 diff)
• _mainLinearBlock.tol (added)
• _outputMissingBlock.tol (modified) (1 diff)
• _stdLinearBlock.tol (modified) (1 diff)

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

@@ -32 +32 @@
   Text _.autodoc.description="Gibbs sampler for model "+_.name;
   Text _MID = "[BSR."+_.name+"] ";
-  NameBlock _.lin.blk = BysMcmc::Bsr::Gibbs::StdLinearBlock
+  NameBlock _.mln.blk = BysMcmc::Bsr::Gibbs::MainLinearBlock
   (
     "BSR.MainLinBlock."+_.name,
@@ -46 +46 @@
     _.modelDef->NoiseDistrib,
     2,
-    _.lin.blk::_.blockDef::_.lastCol+1,
+    _.mln.blk::_.blockDef::_.lastCol+1,
     _config
   );
@@ -79 +79 @@
   Set _.blocks =
   [[
-    _.lin.blk,
+    _.mln.blk,
     _.sig.blk,
     _.arm.blk,
@@ -199 +199 @@
   ////////////////////////////////////////////////////////////////////////////
   {
-  //WriteLn("TRACE"+_MID+" _.lin.blk.init 1");
+  //WriteLn("TRACE"+_MID+" init 1");
     Matrix aux=SetCol(EvalSet(_.modelDef->LinearBlock,Real(BSR.LinearBlock b)
     { b->InitialValue }));
     VMatrix _.YplusV := _.outMis.blk::filter(_.Y);
-  //WriteLn("TRACE"+_MID+" _.lin.blk.init 6");
+  //WriteLn("TRACE"+_MID+" init 6");
     VMatrix _.XplusU := _.inpMis.blk::filter(_.modelDef->X);
-  //WriteLn("TRACE"+_MID+" _.lin.blk.init 2");
+  //WriteLn("TRACE"+_MID+" init 2");
     VMatrix beta0 = findMatchingBeta(Mat2VMat(aux));
-  //WriteLn("TRACE"+_MID+" _.lin.blk.init 3");
-    Real ok = _.lin.blk::initialize(beta0, _.modelDef->A, _.modelDef->a);
-  //WriteLn("TRACE"+_MID+" _.lin.blk.init 4");
-    VMatrix _.beta   := _.lin.blk::_.blockDef::_.store;
+  //WriteLn("TRACE"+_MID+" init 3");
+    Real ok = _.mln.blk::initialize(beta0, _.modelDef->A, _.modelDef->a);
+  //WriteLn("TRACE"+_MID+" init 4");
+    VMatrix _.beta   := _.mln.blk::getStore(0);
     Real DoExternalFilter(0);
     VMatrix sig_0 = _.sig.blk::draw(0, _.Y - _.modelDef->X*_.beta);
-  //WriteLn("TRACE"+_MID+" _.lin.blk.init 7");
+  //WriteLn("TRACE"+_MID+" init 7");
     VMatrix _.filter := _.XplusU*_.beta;
     VMatrix _.noise := _.YplusV-_.filter;
@@ -263 +263 @@
     VMatrix X01 = _.arm.blk::filter(si*_.XplusU); 
   //WriteLn("TRACE"+_MID+" draw.sig.lin 8");
-    VMatrix _.beta := If(Not(_.lin.blk::_.blockDef::_.numCol),Rand(0,1,0,0),
-    {
-      Real _.lin.blk::solve(numSim, Y01, X01);
-      _.lin.blk::draw(numSim, Y01, X01)
-    });
+    VMatrix _.beta := If(Not(_.mln.blk::_.N),Rand(0,1,0,0),
+      _.mln.blk::draw(numSim, Y01, X01));
   //WriteLn("TRACE"+_MID+" draw.sig.lin 9");
     _.beta << sig << arm << im << om
@@ -309 +306 @@
     VMatrix X01 = _.arm.blk::filter(si*_.XplusU); 
   //WriteLn("TRACE"+_MID+" draw.sig.lin 8");
-    VMatrix _.beta := If(Not(_.lin.blk::_.blockDef::_.numCol),Rand(0,1,0,0),
-      _.lin.blk::maxLikelihood(numSim, Y01, X01));
+    VMatrix _.beta := If(Not(_.mln.blk::_.N),Rand(0,1,0,0),
+      _.mln.blk::maxLikelihood(numSim, Y01, X01));
   //WriteLn("TRACE"+_MID+" draw.sig.lin 9");
     _.beta << sig << arm << im << om
@@ -397 +394 @@
     VMatrix _.XplusU := _.inpMis.blk::filter(_.modelDef->X);
   //WriteLn("TRACE"+_MID+" eval 4");
-    VMatrix _.beta := _.lin.blk::_.blockDef::_.store;
+    VMatrix _.beta := _.mln.blk::getStore(0);
     VMatrix fullFilter = _.modelDef->X*_.beta;
     VMatrix fullFilter.full = _.XplusU*_.beta;
@@ -767 +764 @@
   ////////////////////////////////////////////////////////////////////////////
   {
-    NameBlock lin = _.lin.blk::_.blockDef;
+    NameBlock lin = _.mln.blk::_.lin.blk::_.blockDef;
     NameBlock sig = _.sig.blk::_.blockDef;
 

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

@@ -15 +15 @@
   Set _linBlk = Include("_stdLinearBlock.tol");
   Code StdLinearBlock;
+
+  //Linear blocks related methods
+  Set _mnlBlk = Include("_mainLinearBlock.tol");
+  Code MainLinearBlock;
 
   //Sigma block related methods

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

@@ -211 +211 @@
       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);
+      VMatrix _.blockDef::_.store := _.lin.blk::draw(numSim);
       _.blockDef::_.store
     })

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

@@ -192 +192 @@
       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);
+      VMatrix _.blockDef::_.store := _.lin.blk::draw(numSim);
       _.blockDef::_.store
     })

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

@@ -288 +288 @@
   Text _.autodoc.member.draw =
   "Draws a simulation of a linear block in current case, constrained or not";
-  VMatrix draw(Real numSim, VMatrix Y, VMatrix X)
+  VMatrix draw(Real numSim)
   ////////////////////////////////////////////////////////////////////////////
   {
     If(_.blockDef::enabled,
     {
-    //Real solve(numSim,Y,X);
       VMatrix aux = _drawer(numSim);
     //WriteLn("TRACE "+_MID+" draw");