Changeset 6049f18 in tol

Timestamp:

May 14, 2010, 1:14:26 PM (16 years ago)

Author:

Víctor de Buen <vdebuen@…>

Branches:

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

Children:

4f775f9

Parents:

09cf872

Message:

Usando la nomenclatura de Minka
Adaptación para Multiple Try Metropolis

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

File:

• tol_pkg/BysVecLinReg/Multivariate.tol (modified) (9 diffs)

tol_pkg/BysVecLinReg/Multivariate.tol

@@ -62 +62 @@
   //Stacked columns of A (D*Mx1)
   VMatrix _.vec.A = Rand(0,0,0,0); 
-  VMatrix _.vec.A.prior.nu = Rand(0,0,0,0); 
-  VMatrix _.vec.A.prior.S = Rand(0,0,0,0); 
+  VMatrix _.M = Rand(0,0,0,0); 
+  VMatrix _.K = Rand(0,0,0,0); 
   //Average of posterior distribution of A for alpha=0: Y*X'*(X*Xt)^-1 (dxm)
   VMatrix _.nu = Rand(0,0,0,0); 
@@ -69 +69 @@
   VMatrix _.vec.nu = Rand(0,0,0,0); 
   //Y*Y'
-  VMatrix _.YYt = Rand(0,0,0,0); 
+  VMatrix _.Syy = Rand(0,0,0,0); 
+  VMatrix _.Syx = Rand(0,0,0,0); 
+  VMatrix _.Sxx = Rand(0,0,0,0); 
+  VMatrix _.Sy.x = Rand(0,0,0,0); 
+
   //Y*X'*(X*Xt)^-1*X*Y'
   VMatrix _.nuXYt = Rand(0,0,0,0); 
@@ -181 +185 @@
     1==1, 
     {
+      VMatrix _.M := Constant(_.d,_.m,0);
+      VMatrix _.K := Constant(_.m,_.m,0);
       Real change.X(X);
       Real change.Y(Y);
@@ -187 +193 @@
       VMatrix _.A := BysVecLinReg::vec2vmatrix(_.vec.A,Constant(_.d,_.m,1));
       Real change.V(Eye(_.d));
-      VMatrix _.vec.A.prior.nu := _.vec.nu*0;
-      VMatrix _.vec.A.prior.S := _.V*0;
       True
     })
   };
+
 
   ////////////////////////////////////////////////////////////////////////////
@@ -231 +236 @@
     {
       VMatrix _.X := Copy(X); 
-      VMatrix _.X.Lf:=CholeskiFactor(_.X,"XXt");
+      VMatrix _.Sxx := MMtSqr(_.X)+_.K;
+      VMatrix _.X.Lf:=CholeskiFactor(_.Sxx,"X");
       VMatrix _.X.Le := CholeskiSolve(_.X.Lf,
         Convert(_.X.Lf, "Cholmod.R.Sparse"),"P");
@@ -248 +254 @@
     {
       VMatrix _.Y := Copy(Y); 
-      VMatrix _.YYt := _.Y*Tra(_.Y);
+      VMatrix _.Syy := _.Y*Tra(_.Y) + _.M*_.K*Tra(_.M);
+      VMatrix _.Syx := _.Y*Tra(_.X) + _.M*_.K;
       VMatrix XYt =_.X*Tra(_.Y);
-      VMatrix _.nu:=Tra(CholeskiSolve(_.X.Lf,XYt,"PtLLtP"));
-      VMatrix _.nuXYt := _.nu*XYt;
+      VMatrix _.nu:=Tra(CholeskiSolve(_.X.Lf,Tra(_.Syx),"PtLLtP"));
+      VMatrix _.Sy.x := _.Syy - _.nu*Tra(_.Syx);
       VMatrix _.vec.nu:=BysVecLinReg::vmatrix2vec(_.nu);
       True
@@ -317 +324 @@
   ////////////////////////////////////////////////////////////////////////////
   {
-    VMatrix _.E := _.Y - _.A * _.X;
-    //VMatrix S2 = MMtSqr(_.E);
-        VMatrix S2 = _.YYt -_.nuXYt;
-    VMatrix W = S2 + _.W.V0*_.W.N0;
+  //VMatrix S2 = MMtSqr(_.E);
+    VMatrix W = _.Sy.x + _.W.V0*_.W.N0;
     Real freDeg = _.N+_.W.N0;
     VMatrix V = Mat2VMat(RandIWishart(freDeg, VMat2Mat(W)));
@@ -336 +341 @@
 
   ////////////////////////////////////////////////////////////////////////////
-  VMatrix gibbs.draw.vec.A.std(Real numSim)
-  ////////////////////////////////////////////////////////////////////////////
-  {
+  VMatrix gibbs.draw.vec.A.std(Real numSim, VMatrix C01, VMatrix c01)
+  ////////////////////////////////////////////////////////////////////////////
+  {
+  //WriteLn("TRACE gibbs.draw.vec.A.std");
     If(!_.r.C, 
       Gaussian(_.d*_.m,1,0,1),
-      BysVecLinReg::eqConstrainedStdGaussian(_.C01,_.c01))
-  };
-
-  ////////////////////////////////////////////////////////////////////////////
-  VMatrix try.draw.vec.A(Real numSim)
-  ////////////////////////////////////////////////////////////////////////////
-  {
+      BysVecLinReg::eqConstrainedStdGaussian(C01,c01))
+  };
+
+  ////////////////////////////////////////////////////////////////////////////
+  VMatrix try.draw.vec.A_mu_s(Real numSim, VMatrix C01, VMatrix c01, VMatrix nu, Real s)
+  ////////////////////////////////////////////////////////////////////////////
+  {
+  //WriteLn("TRACE try.draw.vec.A_mu_s");
     Real _.totalTry := _.totalTry + 1;
-    VMatrix e = gibbs.draw.vec.A.std(numSim);
+    VMatrix e = gibbs.draw.vec.A.std(numSim,C01,c01)*s;
     VMatrix e.1 = _.VxIm.Le*e;
     VMatrix e.2 = Tra(Tra(e.1)*_.IdxXLi);
-    (_.vec.nu + e.2)/(1+_.alpha)
+    (nu + e.2)/(1+_.alpha)
+  };
+
+  ////////////////////////////////////////////////////////////////////////////
+  VMatrix try.draw.vec.A(Real numSim)
+  ////////////////////////////////////////////////////////////////////////////
+  {
+     try.draw.vec.A_mu_s(numSim,_.C01,_.c01,_.vec.nu, 1)
   };
 
@@ -411 +425 @@
       VMatrix _.A := BysVecLinReg::vec2vmatrix(_.vec.A,_.A)
     });
+    VMatrix _.E := _.Y - _.A * _.X;
     Copy(_.A)
   };