Changeset dd921c2 in tol

Timestamp:

Apr 1, 2009, 3:51:34 PM (18 years ago)

Author:

Víctor de Buen <vdebuen@…>

Branches:

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

Children:

c07b2db

Parents:

24bdea7

Message:

Removing max-likelihood step
Preparing for Metropolis-Hastings method explained in #664

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

File:

• tol/stdlib/math/stat/models/bayesian/bysMcmc/bsr/gibbs/_stdLinearBlock.tol (modified) (11 diffs)

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

@@ -59 +59 @@
   Set _.solution = Copy(Empty);
 
-  Real _changed.X = True;
+  Real changed.X = True;
 
   //Maximum accepted rounding error for Choleski operations
@@ -67 +67 @@
   //Current maximum of rounding error in solution
   Real _.cholSolMaxRoundErr     = 0;
-
+  //Logarithm of density of last generated vector
+  Real _.betaLogDens = -1/0;
   //Stores the associated drawing method
-  Code _drawer = UnknownCode;
+  Code drawer = UnknownCode;
 
   ////////////////////////////////////////////////////////////////////////////
@@ -89 +90 @@
     VMatrix _.a    := a;
   //WriteLn("TRACE "+_MID+" initialize 6 _.a="<<_.a);
-    Code    _drawer   := If(_.R, draw.constrained, draw.unconstrained);
+    Code    drawer   := If(_.R, draw.constrained, draw.unconstrained);
   //WriteLn("TRACE "+_MID+" initialize 7");
     VMatrix (_.blockDef::_.store) := If(!_.R, beta0,
@@ -137 +138 @@
     If(And(design.is.constant,VRows(_.L)), 0,
     {
-      Real _changed.X := If(Or(VColumns(X)!=VColumns(_old.X),
+      Real changed.X := If(Or(VColumns(X)!=VColumns(_old.X),
                                VRows   (X)!=VRows   (_old.X)), -1,
       {
         1+VMatMax(Abs(X-_old.X))-1
       });
-      If(!_changed.X ,0,
+      If(!changed.X ,0,
       {
         Real density = VNonNullCells(X)/(VRows(X)*VColumns(X));
@@ -160 +161 @@
         },        
         {
-        //WriteLn("TRACE"+_MID+" factorize changed="<<_changed.X );
+        //WriteLn("TRACE"+_MID+" factorize changed="<<changed.X );
           Real numErr = Copy(NError);
           VMatrix _.L := If(density<.13, 
@@ -195 +196 @@
   ////////////////////////////////////////////////////////////////////////////
   //Solves the linear system corresponding to the standarized linear regression
-  Real solve(Real numSim, VMatrix Y, VMatrix X)
-  ////////////////////////////////////////////////////////////////////////////
-  {
-  //WriteLn("TRACE "+_MID+" solve 1 X="<< VMat2Mat(X));
+  Real solve(Real numSim_, VMatrix Y_, VMatrix X_)
+  ////////////////////////////////////////////////////////////////////////////
+  {
+    Real numSim = numSim_;
+    VMatrix Y = Y_;
+    VMatrix X = X_;
+  //WriteLn("TRACE "+_MID+" solve 1 ");
     Real Do factorize(X);
     Real If(VRows(_.L)<=0, CMsg::Coded::show
@@ -235 +239 @@
   Text _.autodoc.member.draw.unconstrained =
   "Method of generation of a linear block in unconstrained case";
-  VMatrix draw.unconstrained(Real numSim)
-  ////////////////////////////////////////////////////////////////////////////
-  {
-  //WriteLn("TRACE "+_MID+" draw.unconstrained 2");
-    VMatrix   b1 = Gaussian(_.N,1,0,1);
-  //WriteLn("TRACE "+_MID+" draw.unconstrained 3 b1="<<Matrix(VMat2Mat(Tra(b1))));
-    VMatrix   b2 = CholeskiSolve(_.L,b1,"LtP");
-  //WriteLn("TRACE "+_MID+" draw.unconstrained 4 _.L="<<Matrix(VMat2Mat(_.L)));
-  //WriteLn("TRACE "+_MID+" draw.unconstrained 4 b2="<<Matrix(VMat2Mat(Tra(b2))));
-    VMatrix  (_.blockDef::_.store) := _.solution::nu+b2;
-  //WriteLn("TRACE "+_MID+" draw.unconstrained 5");
+  VMatrix draw.unconstrained(Real numSim, VMatrix nu)
+  ////////////////////////////////////////////////////////////////////////////
+  {
+    VMatrix z = Gaussian(_.N,1,0,1);
+    Real _.betaLogDens := -0.5*VMatDat(MtMSqr(z),1,1);
+    VMatrix  (_.blockDef::_.store):=nu+CholeskiSolve(_.L,z,"LtP");
     Copy(_.blockDef::_.store)
   };
@@ -252 +251 @@
   Text _.autodoc.member.draw.constrained =
   "Method of generation of a linear block in constrained case";
-  VMatrix draw.constrained(Real numSim)
-  ////////////////////////////////////////////////////////////////////////////
-  {
-    VMatrix nu = _.solution::nu;
+  VMatrix draw.constrained(Real numSim, VMatrix nu)
+  ////////////////////////////////////////////////////////////////////////////
+  {
   //WriteLn("TRACE "+_MID+" draw.constrained 1.1 Max(A*b0-a)="<<VMatMax(_.A*_.blockDef::_.store - a));
     VMatrix d   = _.a - _.A*nu;
-    VMatrix If(_changed.X, _.D  := Tra(CholeskiSolve(_.L,_.At,"PtL")));
+    VMatrix If(changed.X, _.D  := Tra(CholeskiSolve(_.L,_.At,"PtL")));
     VMatrix z0  = Tra(Tra(_.blockDef::_.store-nu)*_.Ls);
   //WriteLn("TRACE "+_MID+" draw.constrained 1.2 Max(D*z0-d)="<<VMatMax(_.D*z0 - d));
@@ -281 +279 @@
       1
     });
+    Real _.betaLogDens := -0.5*VMatDat(MtMSqr(z),1,1);
     VMatrix  (_.blockDef::_.store) := nu + CholeskiSolve(_.L,z,"LtP");
     Copy(_.blockDef::_.store)
@@ -293 +292 @@
     If(_.blockDef::enabled,
     {
-      VMatrix aux = _drawer(numSim);
+      VMatrix aux = drawer(numSim, _.solution::nu);
     //WriteLn("TRACE "+_MID+" draw");
       Real _.blockDef::checkNonFiniteValues(numSim);
@@ -304 +303 @@
 
   ////////////////////////////////////////////////////////////////////////////
-  Text _.autodoc.member.maxLikelihood =
-  "Method of contitioned maximum likelihood";
-  VMatrix maxLikelihood(Real numSim, VMatrix Y, VMatrix X)
-  ////////////////////////////////////////////////////////////////////////////
-  {
-    Real solve(numSim,Y,X);
-    VMatrix aux = If(!_.R, 
-    {
-      _.solution::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)
-      })
-    });
-    VMatrix  (_.blockDef::_.store) := aux;
-  //WriteLn("TRACE "+_MID+" maxLikelihood aux="<<Matrix(Tra(VMat2Mat(aux))));
-    Real _.blockDef::checkNonFiniteValues(numSim);
-    Copy(_.blockDef::_.store)
-  };
-
-
-  ////////////////////////////////////////////////////////////////////////////
   Text _.autodoc.member.setStore =
   "Sets _.store to the specified row vector of values";