Index: tol_pkg/BysVecLinReg/BysVecLinReg.tol
===================================================================
--- tol_pkg/BysVecLinReg/BysVecLinReg.tol	(revision fd471abab5c7f913685f1124178262bea3223422)
+++ tol_pkg/BysVecLinReg/BysVecLinReg.tol	(revision 412906ecc4a1cd831cca41642d3bd0ba84fd37ae)
@@ -9,9 +9,9 @@
 Text _.autodoc.url = "";
 Set _.autodoc.keys = [["MCMC","Markov chain","Bayesian","Metropolis",
-  "Hastings","multivariant","linear","regression","inequality","constraints"]];
+  "Hastings","multivariate","linear","regression","inequality","constraints"]];
 Set _.autodoc.authors = [[
    "vdebuen@tol-project.org"]];
 
-#Embed "Multivariant.tol"
+#Embed "Multivariate.tol"
 
 ]];
Index: l_pkg/BysVecLinReg/Multivariant.tol
===================================================================
--- tol_pkg/BysVecLinReg/Multivariant.tol	(revision fd471abab5c7f913685f1124178262bea3223422)
+++ 	(revision )
@@ -1,286 +1,0 @@
-/* -*- mode: C++ -*- */
-//////////////////////////////////////////////////////////////////////////////
-// FILE   : Multivariant.tol
-// PURPOSE: Defines Class BysVecLinReg::@Multivariant
-//////////////////////////////////////////////////////////////////////////////
-
-//////////////////////////////////////////////////////////////////////////////
-//Multivariant linear regression
-// Y = A*X + E
-//where
-//   Y is the known output matrix DxN
-//   X is the known input matrix MxN
-//   A is the unknown cofficient matrix NxD
-//   E is the unknown residuals matrix DxN with distribution
-//     E[.,t] ~ N(0,V) for all t = 1 ... N
-//     V has a prior inverse Wishart W(V0,N0)
-//
-Class @Multivariant
-//////////////////////////////////////////////////////////////////////////////
-{
-  Text _.desc;
-  Set _.outputNames;
-  Set _.inputNames;
-  ////////////////////////////////////////////////////////////////////////////
-  //Definition members
-  ////////////////////////////////////////////////////////////////////////////
-  VMatrix _.Y = Rand(0,0,0,0);    //Multivariant output matrix DxN 
-  VMatrix _.X = Rand(0,0,0,0);    //Multivariant input matrix MxN
-  Real _.alpha = 0;               //Prior invariant scaling factor
-  Real _.W.weight = 0;            //Weight of prior Wishart freedom degrees
-  Real _.W.N0 = 0;                //Prior Wishart freedom degrees
-  VMatrix _.W.V0 = Rand(0,0,0,0); //Prior Wishart scale factor matrix
-  
-  ////////////////////////////////////////////////////////////////////////////
-  //Auxiliar members
-  ////////////////////////////////////////////////////////////////////////////
-  //Data length
-  Real _.N = ?; 
-  //Number of multivariant output vectors
-  Real _.d = ?; 
-  //Number of multivariantinput vectors
-  Real _.m = ?; 
-  //multivariant residuals (dxN)
-  VMatrix _.E = Rand(0,0,0,0);
-  //Instant covariance of multivariant residuals (dxd)
-  VMatrix _.V = Rand(0,0,0,0);
-  //Regression multivariant coefficients (dxm)
-  VMatrix _.A = Rand(0,0,0,0); 
-  //Stacked columns of A (D*Mx1)
-  VMatrix _.vec.A = 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); 
-  //Stacked columns of nu ((d*m)x1)
-  VMatrix _.vec.nu = Rand(0,0,0,0); 
-  //Y*Y'
-  VMatrix _.YYt = Rand(0,0,0,0); 
-  //Y*X'*(X*Xt)^-1*X*Y'
-  VMatrix _.nuXYt = Rand(0,0,0,0); 
-  //Cholesky decomposition of X*X' as sparse factor (mxm)
-  VMatrix _.X.Lf = Rand(0,0,0,0);
-  //Cholesky decomposition of X*X' as sparse explicit (mxm)
-  VMatrix _.X.Le = Rand(0,0,0,0);
-  //Inverse of Cholesky decomposition of X*X' as sparse explicit (mxm)
-  VMatrix _.X.Li = Rand(0,0,0,0);
-  //Kronecker product of identity of order d by X.L ((d*m)x(d*m))
-  VMatrix _.IdxXL = Rand(0,0,0,0);
-  //Kronecker product of identity of order d by X.Li ((d*m)x(d*m))
-  VMatrix _.IdxXLi = Rand(0,0,0,0);
-  //Kronecker product of V by identity of order m ((d*m)x(d*m))
-  VMatrix _.VxIm = Rand(0,0,0,0);
-  //Choleski decomposition of VxIm as sparse factor((d*m)x(d*m))
-  VMatrix _.VxIm.Lf = Rand(0,0,0,0);
-  //Choleski decomposition of VxIm as sparse explicit((d*m)x(d*m))
-  VMatrix _.VxIm.Le = Rand(0,0,0,0);
- 
-  ////////////////////////////////////////////////////////////////////////////
-  Static VMatrix vmatrix2vec(VMatrix M)
-  ////////////////////////////////////////////////////////////////////////////
-  {
-    Matrix Mt = VMat2Mat(M,True);
-    Mat2VMat(Tra(GetNumeric([[ Mt ]])))
-  };
-
-  ////////////////////////////////////////////////////////////////////////////
-  Static VMatrix vec2vmatrix(VMatrix vec, VMatrix M)
-  ////////////////////////////////////////////////////////////////////////////
-  {
-    Matrix Mt = VMat2Mat(M,True);
-    Mat2VMat(DeepCopy([[Mt]],VMat2Mat(vec))[1],True)
-  };
-
-  ////////////////////////////////////////////////////////////////////////////
-  Real check.dim(Text method, VMatrix M, Real nRow, Real nCol)
-  ////////////////////////////////////////////////////////////////////////////
-  {
-    If(And(VRows(M)==nRow, VColumns(M)==nCol), True,
-    {
-      WriteLn("[BysVecLinReg::@Multivariant::"+method+"] "
-              "Incongruent dimensions of _.W.V0 that should be "
-              "("<<nRow+","<<nCol+") instead of "
-              "("<<VRows(M)+","<<VColumns(M)+")","E"); 
-      False
-    })
-  };
-
-  ////////////////////////////////////////////////////////////////////////////
-  Real set.outputInput(VMatrix Y, VMatrix X)
-  ////////////////////////////////////////////////////////////////////////////
-  {
-    Real _.d := VRows(Y);
-    Real _.m := VRows(X);
-    Real _.N := VColumns(X);
-    Case(
-    VColumns(Y)!=_.N, 
-    {
-      WriteLn("[BysVecLinReg::@Multivariant::set.outputInput] "
-              "Incongruent rows data length between input "
-              "("<<VRows(X)+","<<VColumns(X)+") "
-              "and output "
-              "("<<VRows(Y)+","<<VColumns(Y)+")","E"); 
-      False
-    },
-    _.m > _.N, 
-    {
-      WriteLn("[BysVecLinReg::@Multivariant::set.outputInput] "
-              "Too much inputs for data length "<<_.N+"<"<<_.m,"E"); 
-      False
-    },
-    1==1, 
-    {
-      change.X(X);
-      change.Y(Y);
-      change.alpha(0);
-      VMatrix _.vec.A := (_.vec.nu)/(1+_.alpha);
-      VMatrix _.A := @Multivariant::vec2vmatrix(_.vec.A,Constant(_.d,_.m,1));
-      change.V(Eye(_.d));
-      True
-    })
-  };
-
-  ////////////////////////////////////////////////////////////////////////////
-  Static @Multivariant New(  
-    Text desc,
-    Set outputNames,
-    Set inputNames,
-    VMatrix Y, 
-    VMatrix X)
-  ////////////////////////////////////////////////////////////////////////////
-  {
-    @Multivariant new = [[
-      Text _.desc = desc;
-      Set _.outputNames = outputNames;
-      Set _.inputNames = inputNames]];
-    Real new::set.outputInput(Y,X);
-    new
-  };
-
-  ////////////////////////////////////////////////////////////////////////////
-  Real set.prior.V.inverseWishart(Real W.weight, VMatrix W.V0)
-  ////////////////////////////////////////////////////////////////////////////
-  {
-    If(!check.dim("set.prior.V.inverseWishart",W.V0,_.d,_.d),False,
-    {
-      Real _.W.N0 := Case(
-      W.weight==0,0,
-      W.weight>=0.99,_.N*99,
-      1==1,
-      {
-        Max(1, _.N * W.weight / (1-W.weight))
-      });
-      Real _.W.weight := _.W.N0 / (_.W.N0+_.N);
-      VMatrix _.W.V0 := Copy(W.V0);
-      True
-    })
-  };
-
-  ////////////////////////////////////////////////////////////////////////////
-  Real change.X(VMatrix X)
-  ////////////////////////////////////////////////////////////////////////////
-  {
-    If(!check.dim("change.X",X,_.m,_.N),False,
-    {
-      VMatrix _.X := Copy(X); 
-      VMatrix _.X.Lf:=CholeskiFactor(_.X,"XXt");
-      VMatrix _.X.Le := CholeskiSolve(_.X.Lf,
-        Convert(_.X.Lf, "Cholmod.R.Sparse"),"P");
-      VMatrix _.X.Li := CholeskiSolve(_.X.Lf,Eye(_.m),"PtL");
-      VMatrix _.IdxXL := KroneckerProduct(Eye(_.d),_.X.Le);
-      VMatrix _.IdxXLi := KroneckerProduct(Eye(_.d),_.X.Li);
-      True
-    })
-  };
-
-  ////////////////////////////////////////////////////////////////////////////
-  Real change.Y(VMatrix Y)
-  ////////////////////////////////////////////////////////////////////////////
-  {
-    If(!check.dim("change.Y",Y,_.d,_.N),False,
-    {
-      VMatrix _.Y := Copy(Y); 
-      VMatrix _.YYt := MMtSqr(_.Y);
-      VMatrix XYt =_.X*Tra(_.Y);
-      VMatrix _.nu:=Tra(CholeskiSolve(_.X.Lf,XYt,"PtLLtP"));
-      VMatrix _.nuXYt := _.nu*XYt;
-      VMatrix _.vec.nu:=@Multivariant::vmatrix2vec(_.nu);
-      True
-    })
-  };
-
-  ////////////////////////////////////////////////////////////////////////////
-  Real change.alpha(Real alpha)
-  ////////////////////////////////////////////////////////////////////////////
-  {
-    If(alpha<0,
-    {
-      WriteLn("[BysVecLinReg::@Multivariant::change.alpha] "
-              "scaling invariant parameter alpha must be non negative","E"); 
-      False
-    },
-    {
-      Real _.alpha := alpha;
-      True
-    })
-  };
-
-  ////////////////////////////////////////////////////////////////////////////
-  Real change.V(VMatrix V)
-  ////////////////////////////////////////////////////////////////////////////
-  {
-    If(!check.dim("change.V",V,_.d,_.d),False,
-    {
-      VMatrix _.V := Copy(V); 
-      VMatrix _.VxIm := KroneckerProduct(_.V, Eye(_.m));
-      VMatrix _.VxIm.Lf := CholeskiFactor(_.VxIm,"X");
-      VMatrix _.VxIm.Le := CholeskiSolve(_.VxIm.Lf,
-        Convert(_.VxIm.Lf, "Cholmod.R.Sparse"),"P");
-      VMatrix _.VxIm.Lf := CholeskiFactor(_.VxIm,"X");
-      VMatrix _.VxIm.Le := CholeskiSolve(_.VxIm.Lf,
-        Convert(_.VxIm.Lf, "Cholmod.R.Sparse"),"P");
-      True
-    })
-  };
-
-  ////////////////////////////////////////////////////////////////////////////
-  Real set.maxEvidence.alpha(Real void)
-  ////////////////////////////////////////////////////////////////////////////
-  {
-    VMatrix V.L = CholeskiFactor(_.V,"X");
-    VMatrix aux = CholeskiSolve(V.L,_.nuXYt,"PtLLtP");
-    Real trace = VMatSum(SubDiag(aux));
-    Real md = _.m*_.d;
-    change.alpha(md/(trace-md))
-  };
-
-  ////////////////////////////////////////////////////////////////////////////
-  VMatrix gibbs.draw.V(Real numSim)
-  ////////////////////////////////////////////////////////////////////////////
-  {
-    VMatrix E = _.Y - _.A * _.X;
-    VMatrix S2 = MMtSqr(E);
-    Real change.V(Mat2VMat(RandIWishart(_.N+_.W.N0, VMat2Mat(S2+_.W.V0*_.W.N0))));
-    Copy(_.V)
-  };
-
-  ////////////////////////////////////////////////////////////////////////////
-  VMatrix gibbs.draw.vec.A(Real numSim)
-  ////////////////////////////////////////////////////////////////////////////
-  {
-    VMatrix e = Gaussian(_.d*_.m,1,0,1);
-    VMatrix e.1 = _.VxIm.Le*e;
-    VMatrix e.2 = Tra(Tra(e.1)*_.IdxXLi);
-    VMatrix _.vec.A := (_.vec.nu + e.2)/(1+_.alpha);
-    _.vec.A
-  };
-
-  ////////////////////////////////////////////////////////////////////////////
-  VMatrix gibbs.draw.A(Real numSim)
-  ////////////////////////////////////////////////////////////////////////////
-  {
-    VMatrix gibbs.draw.vec.A(numSim);
-    VMatrix _.A := @Multivariant::vec2vmatrix(_.vec.A,_.A);
-    Copy(_.A)
-  }
-
-}
-
Index: tol_pkg/BysVecLinReg/Multivariate.tol
===================================================================
--- tol_pkg/BysVecLinReg/Multivariate.tol	(revision 412906ecc4a1cd831cca41642d3bd0ba84fd37ae)
+++ tol_pkg/BysVecLinReg/Multivariate.tol	(revision 412906ecc4a1cd831cca41642d3bd0ba84fd37ae)
@@ -0,0 +1,289 @@
+/* -*- mode: C++ -*- */
+//////////////////////////////////////////////////////////////////////////////
+// FILE   : Multivariate.tol
+// PURPOSE: Defines Class BysVecLinReg::@Multivariate
+//////////////////////////////////////////////////////////////////////////////
+
+//////////////////////////////////////////////////////////////////////////////
+//Multivariate linear regression
+// Y = A*X + E
+//where
+//   Y is the known output matrix DxN
+//   X is the known input matrix MxN
+//   A is the unknown cofficient matrix NxD
+//   E is the unknown residuals matrix DxN with distribution
+//     E[.,t] ~ N(0,V) for all t = 1 ... N
+//     V has a prior inverse Wishart W(V0,N0)
+//
+Class @Multivariate
+//////////////////////////////////////////////////////////////////////////////
+{
+  Text _.desc;
+  Set _.outputNames;
+  Set _.inputNames;
+  ////////////////////////////////////////////////////////////////////////////
+  //Definition members
+  ////////////////////////////////////////////////////////////////////////////
+  VMatrix _.Y = Rand(0,0,0,0);    //Multivariate output matrix DxN 
+  VMatrix _.X = Rand(0,0,0,0);    //Multivariate input matrix MxN
+  Real _.alpha = 0;               //Prior invariant scaling factor
+  Real _.W.weight = 0;            //Weight of prior Wishart freedom degrees
+  Real _.W.N0 = 0;                //Prior Wishart freedom degrees
+  VMatrix _.W.V0 = Rand(0,0,0,0); //Prior Wishart scale factor matrix
+  
+  ////////////////////////////////////////////////////////////////////////////
+  //Auxiliar members
+  ////////////////////////////////////////////////////////////////////////////
+  //Data length
+  Real _.N = ?; 
+  //Number of multivariate output vectors
+  Real _.d = ?; 
+  //Number of multivariantinput vectors
+  Real _.m = ?; 
+  //multivariate residuals (dxN)
+  VMatrix _.E = Rand(0,0,0,0);
+  //Instant covariance of multivariate residuals (dxd)
+  VMatrix _.V = Rand(0,0,0,0);
+  //Regression multivariate coefficients (dxm)
+  VMatrix _.A = Rand(0,0,0,0); 
+  //Stacked columns of A (D*Mx1)
+  VMatrix _.vec.A = 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); 
+  //Stacked columns of nu ((d*m)x1)
+  VMatrix _.vec.nu = Rand(0,0,0,0); 
+  //Y*Y'
+  VMatrix _.YYt = Rand(0,0,0,0); 
+  //Y*X'*(X*Xt)^-1*X*Y'
+  VMatrix _.nuXYt = Rand(0,0,0,0); 
+  //Cholesky decomposition of X*X' as sparse factor (mxm)
+  VMatrix _.X.Lf = Rand(0,0,0,0);
+  //Cholesky decomposition of X*X' as sparse explicit (mxm)
+  VMatrix _.X.Le = Rand(0,0,0,0);
+  //Inverse of Cholesky decomposition of X*X' as sparse explicit (mxm)
+  VMatrix _.X.Li = Rand(0,0,0,0);
+  //Kronecker product of identity of order d by X.L ((d*m)x(d*m))
+  VMatrix _.IdxXL = Rand(0,0,0,0);
+  //Kronecker product of identity of order d by X.Li ((d*m)x(d*m))
+  VMatrix _.IdxXLi = Rand(0,0,0,0);
+  //Kronecker product of V by identity of order m ((d*m)x(d*m))
+  VMatrix _.VxIm = Rand(0,0,0,0);
+  //Choleski decomposition of VxIm as sparse factor((d*m)x(d*m))
+  VMatrix _.VxIm.Lf = Rand(0,0,0,0);
+  //Choleski decomposition of VxIm as sparse explicit((d*m)x(d*m))
+  VMatrix _.VxIm.Le = Rand(0,0,0,0);
+ 
+  ////////////////////////////////////////////////////////////////////////////
+  Static VMatrix vmatrix2vec(VMatrix M)
+  ////////////////////////////////////////////////////////////////////////////
+  {
+    Matrix Mt = VMat2Mat(M,True);
+    Mat2VMat(Tra(GetNumeric([[ Mt ]])))
+  };
+
+  ////////////////////////////////////////////////////////////////////////////
+  Static VMatrix vec2vmatrix(VMatrix vec, VMatrix M)
+  ////////////////////////////////////////////////////////////////////////////
+  {
+    Matrix Mt = VMat2Mat(M,True);
+    Mat2VMat(DeepCopy([[Mt]],VMat2Mat(vec))[1],True)
+  };
+
+  ////////////////////////////////////////////////////////////////////////////
+  Real check.dim(Text method, VMatrix M, Real nRow, Real nCol)
+  ////////////////////////////////////////////////////////////////////////////
+  {
+    If(And(VRows(M)==nRow, VColumns(M)==nCol), True,
+    {
+      WriteLn("[BysVecLinReg::@Multivariate::"+method+"] "
+              "Incongruent dimensions of _.W.V0 that should be "
+              "("<<nRow+","<<nCol+") instead of "
+              "("<<VRows(M)+","<<VColumns(M)+")","E"); 
+      False
+    })
+  };
+
+  ////////////////////////////////////////////////////////////////////////////
+  Real set.outputInput(VMatrix Y, VMatrix X)
+  ////////////////////////////////////////////////////////////////////////////
+  {
+    Real _.d := VRows(Y);
+    Real _.m := VRows(X);
+    Real _.N := VColumns(X);
+    Case(
+    VColumns(Y)!=_.N, 
+    {
+      WriteLn("[BysVecLinReg::@Multivariate::set.outputInput] "
+              "Incongruent rows data length between input "
+              "("<<VRows(X)+","<<VColumns(X)+") "
+              "and output "
+              "("<<VRows(Y)+","<<VColumns(Y)+")","E"); 
+      False
+    },
+    _.m > _.N, 
+    {
+      WriteLn("[BysVecLinReg::@Multivariate::set.outputInput] "
+              "Too much inputs for data length "<<_.N+"<"<<_.m,"E"); 
+      False
+    },
+    1==1, 
+    {
+      change.X(X);
+      change.Y(Y);
+      change.alpha(0);
+      VMatrix _.vec.A := (_.vec.nu)/(1+_.alpha);
+      VMatrix _.A := @Multivariate::vec2vmatrix(_.vec.A,Constant(_.d,_.m,1));
+      change.V(Eye(_.d));
+      True
+    })
+  };
+
+  ////////////////////////////////////////////////////////////////////////////
+  Real set.prior.V.none(Real void)
+  ////////////////////////////////////////////////////////////////////////////
+  {
+    Real _.W.N0 := 0;
+    Real _.W.weight := 0;
+    VMatrix _.W.V0 := Constant(_.d,_.d,0);
+    VMatrix gibbs.draw.V(0);
+    True
+  };
+
+  ////////////////////////////////////////////////////////////////////////////
+  Real set.prior.V.inverseWishart(Real W.weight, VMatrix W.V0)
+  ////////////////////////////////////////////////////////////////////////////
+  {
+    If(!check.dim("set.prior.V.inverseWishart",W.V0,_.d,_.d),False,
+    {
+      Real _.W.N0 := Case(
+      W.weight==0,0,
+      W.weight>=0.99,_.N*99,
+      1==1,
+      {
+        Max(1, _.N * W.weight / (1-W.weight))
+      });
+      Real _.W.weight := _.W.N0 / (_.W.N0+_.N);
+      VMatrix _.W.V0 := Copy(W.V0);
+      VMatrix gibbs.draw.V(0);
+      True
+    })
+  };
+
+  ////////////////////////////////////////////////////////////////////////////
+  Real change.X(VMatrix X)
+  ////////////////////////////////////////////////////////////////////////////
+  {
+    If(!check.dim("change.X",X,_.m,_.N),False,
+    {
+      VMatrix _.X := Copy(X); 
+      VMatrix _.X.Lf:=CholeskiFactor(_.X,"XXt");
+      VMatrix _.X.Le := CholeskiSolve(_.X.Lf,
+        Convert(_.X.Lf, "Cholmod.R.Sparse"),"P");
+      VMatrix _.X.Li := CholeskiSolve(_.X.Lf,Eye(_.m),"PtL");
+      VMatrix _.IdxXL := KroneckerProduct(Eye(_.d),_.X.Le);
+      VMatrix _.IdxXLi := KroneckerProduct(Eye(_.d),_.X.Li);
+      True
+    })
+  };
+
+  ////////////////////////////////////////////////////////////////////////////
+  Real change.Y(VMatrix Y)
+  ////////////////////////////////////////////////////////////////////////////
+  {
+    If(!check.dim("change.Y",Y,_.d,_.N),False,
+    {
+      VMatrix _.Y := Copy(Y); 
+      VMatrix _.YYt := MMtSqr(_.Y);
+      VMatrix XYt =_.X*Tra(_.Y);
+      VMatrix _.nu:=Tra(CholeskiSolve(_.X.Lf,XYt,"PtLLtP"));
+      VMatrix _.nuXYt := _.nu*XYt;
+      VMatrix _.vec.nu:=@Multivariate::vmatrix2vec(_.nu);
+      True
+    })
+  };
+
+  ////////////////////////////////////////////////////////////////////////////
+  Real change.alpha(Real alpha)
+  ////////////////////////////////////////////////////////////////////////////
+  {
+    If(alpha<0,
+    {
+      WriteLn("[BysVecLinReg::@Multivariate::change.alpha] "
+              "scaling invariant parameter alpha must be non negative","E"); 
+      False
+    },
+    {
+      Real _.alpha := alpha;
+      True
+    })
+  };
+
+  ////////////////////////////////////////////////////////////////////////////
+  Real change.V(VMatrix V)
+  ////////////////////////////////////////////////////////////////////////////
+  {
+    If(!check.dim("change.V",V,_.d,_.d),False,
+    {
+      VMatrix _.V := (V+Tra(V))/2; 
+      VMatrix _.VxIm := KroneckerProduct(_.V, Eye(_.m));
+      VMatrix _.VxIm.Lf := CholeskiFactor(_.VxIm,"X");
+      VMatrix _.VxIm.Le := CholeskiSolve(_.VxIm.Lf,
+        Convert(_.VxIm.Lf, "Cholmod.R.Sparse"),"P");
+      VMatrix _.VxIm.Lf := CholeskiFactor(_.VxIm,"X");
+      VMatrix _.VxIm.Le := CholeskiSolve(_.VxIm.Lf,
+        Convert(_.VxIm.Lf, "Cholmod.R.Sparse"),"P");
+      True
+    })
+  };
+
+  ////////////////////////////////////////////////////////////////////////////
+  Real set.maxEvidence.alpha(Real void)
+  ////////////////////////////////////////////////////////////////////////////
+  {
+    VMatrix V.L = CholeskiFactor(_.V,"X");
+    VMatrix aux = CholeskiSolve(V.L,_.nuXYt,"PtLLtP");
+    Real trace = VMatSum(SubDiag(aux));
+    Real md = _.m*_.d;
+    change.alpha(md/(trace-md))
+  };
+
+  ////////////////////////////////////////////////////////////////////////////
+  VMatrix gibbs.draw.V(Real numSim)
+  ////////////////////////////////////////////////////////////////////////////
+  {
+    VMatrix _.E := _.Y - _.A * _.X;
+    VMatrix S2 = MMtSqr(_.E);
+    Real change.V(Mat2VMat(RandIWishart(_.N+_.W.N0, VMat2Mat(S2+_.W.V0*_.W.N0))));
+    Copy(_.V)
+  };
+
+  ////////////////////////////////////////////////////////////////////////////
+  VMatrix gibbs.draw.vec.V(Real numSim)
+  ////////////////////////////////////////////////////////////////////////////
+  {
+    VMatrix gibbs.draw.V(numSim);
+    @Multivariate::vmatrix2vec(_.V)
+  };
+
+  ////////////////////////////////////////////////////////////////////////////
+  VMatrix gibbs.draw.vec.A(Real numSim)
+  ////////////////////////////////////////////////////////////////////////////
+  {
+    VMatrix e = Gaussian(_.d*_.m,1,0,1);
+    VMatrix e.1 = _.VxIm.Le*e;
+    VMatrix e.2 = Tra(Tra(e.1)*_.IdxXLi);
+    VMatrix _.vec.A := (_.vec.nu + e.2)/(1+_.alpha);
+    _.vec.A
+  };
+
+  ////////////////////////////////////////////////////////////////////////////
+  VMatrix gibbs.draw.A(Real numSim)
+  ////////////////////////////////////////////////////////////////////////////
+  {
+    VMatrix gibbs.draw.vec.A(numSim);
+    VMatrix _.A := @Multivariate::vec2vmatrix(_.vec.A,_.A);
+    Copy(_.A)
+  }
+
+}
+
Index: tol_pkg/BysVecLinReg/test/test_0002/test.tol
===================================================================
--- tol_pkg/BysVecLinReg/test/test_0002/test.tol	(revision fd471abab5c7f913685f1124178262bea3223422)
+++ tol_pkg/BysVecLinReg/test/test_0002/test.tol	(revision 412906ecc4a1cd831cca41642d3bd0ba84fd37ae)
@@ -8,5 +8,5 @@
 Text link_ = "";
 Text summary_ = 
-"test_0001 multivariate linear model with unknown covariance matrix" ;
+"test_0002 multivariate linear model with unknown covariance matrix" ;
 
 Real numErr0 = Copy(NError);