Ticket #507: test3.tol

//////////////////////////////////////////////////////////////////////////////
// FILE   : test.tol
// PURPOSE: 
//////////////////////////////////////////////////////////////////////////////

//////////////////////////////////////////////////////////////////////////////
// INCLUDES
//////////////////////////////////////////////////////////////////////////////
//Set Include(PathIniSadd(0));

//////////////////////////////////////////////////////////////////////////////
// FUNCTIONS
//////////////////////////////////////////////////////////////////////////////

Matrix Diag2Col(Matrix M)
{
  Matrix SetCol(Set For(1,Real Min(Rows(M),Columns(M)), Real (Real r)
      {
        Real MatDat(M,r,r)
      }))
};

Set CheckConstraints(Matrix sample, Matrix B, Matrix rhs)
{
  Set CC = For(1, Columns(sample), Matrix(Real c) {SubCol(sample,[[c]])});
  Set Select(CC, Real(Matrix beta) {
    Matrix check = LE(B*beta,rhs);
    Real NE(Rows(check), MatSum(check))
  })
};

//////////////////////////////////////////////////////////////////////////////
// PROCEDURES
//////////////////////////////////////////////////////////////////////////////

//Set info  = waredrobe[1];
//Real Ois.Store(info, "test3.oza"); 
Set info2 = Ois.Load("test3.oza")[1];


Matrix nu2 = info2["nu"];
Matrix S2 = info2["S"];
Matrix BroNew2 = info2["BroNew"];
Matrix broNew2 = info2["broNew"];
Matrix paramIni2 = info2["paramIni"];

Set check_beta0 = CheckConstraints(paramIni2, BroNew2, broNew2);

Text If(Card(check_beta0),
  WriteLn("beta0 no cumple"),WriteLn("beta0 ok"));

Real Size = 3;
Real Burning = 1;

Set svd = SVD(S2);
Matrix Choleski(S2);

Matrix D1 = BroNew2*svd[1];
//Set View([[D1]], "Std");

Matrix D = D1*SqRt(svd[2]);

//Set View([[D]], "Std");

Matrix markovChain0 = 
GibbsConstrainedMNormal(nu2, S2,
                        BroNew2, broNew2, 1, Size, Burning, paramIni2);


Set check_mcmc0 = CheckConstraints(Tra(markovChain0), BroNew2, broNew2);

Text If(Card(check_mcmc0),
  WriteLn("markovChain0 no cumple"),WriteLn("markovChain0 ok"));

Matrix markovChain1=
GibbsConstrainedMNormal(nu2, [["FACT_SVD", svd[1], SqRt(Diag2Col(svd[2]))]],
                        BroNew2, broNew2, 1, Size, Burning, paramIni2);

Set check_mcmc1 = CheckConstraints(Tra(markovChain1), BroNew2, broNew2);

Text If(Card(check_mcmc1),
  WriteLn("markovChain1 no cumple"),WriteLn("markovChain1 ok"));

Matrix markovChain2=
GibbsConstrainedMNormal(nu2, [["FULL_SVD", S2]],
                        BroNew2, broNew2, 1, Size, Burning, paramIni2);

Set check_mcmc2 = CheckConstraints(Tra(markovChain2), BroNew2, broNew2);

Text If(Card(check_mcmc2),
  WriteLn("markovChain2 no cumple"),WriteLn("markovChain2 ok"));