Changeset 0c19991 in tol


Ignore:
Timestamp:
Aug 12, 2009, 10:34:58 AM (17 years ago)
Author:
Víctor de Buen <vdebuen@…>
Branches:
1941-docker-image-to-test, devel, master
Children:
d3756e3
Parents:
a2d7c82
Message:

Preparing StdLinearBlock and inherited classes for partial sampling

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

Location:
tol/stdlib/math/stat/models/bayesian/bysMcmc/bsr/gibbs
Files:
4 edited

Legend:

Unmodified
Added
Removed

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

    ra2d7c82 r0c19991  
    2222  If(!Card(inputMissingBlock),
    2323  BysMcmc::Bsr::Gibbs::StdLinearBlock.Default(name, numBlock, firstCol, config),
    24   BysMcmc::Bsr::Gibbs::StdLinearFullBlock blk = [[
     24  BysMcmc::Bsr::Gibbs::StdLinearBlock blk = [[
    2525  //Given initial values of block variables
    2626  VMatrix _.store_ = Mat2VMat(SetMat(Extract(inputMissingBlock, 6)));

  • tol/stdlib/math/stat/models/bayesian/bysMcmc/bsr/gibbs/_mainLinearBlock.tol

    ra2d7c82 r0c19991  
    3030  If(!Card(varNames),
    3131  BysMcmc::Bsr::Gibbs::StdLinearBlock.Default(name,numBlock,firstCol,config),
    32   BysMcmc::Bsr::Gibbs::StdLinearFullBlock blk = [[
     32  BysMcmc::Bsr::Gibbs::StdLinearBlock blk = [[
    3333
    3434  ////////////////////////////////////////////////////////////////////////////
    35   Real _buildWorkSpace(VMatrix  Y,    //Non filtered output
    36                        VMatrix  X,    //Non filtered input
    37                        BysMcmc::BlockSampler arm.blk) //ArimaBlock filter
     35  Real _do_buildWorkSpace(VMatrix Y,    //Non filtered output
     36                          VMatrix X,    //Non filtered input
     37                          BysMcmc::BlockSampler arm.blk) //ArimaBlock filter
    3838  ////////////////////////////////////////////////////////////////////////////
    3939  {
     40  //WriteLn("TRACE"+_MID+" MainLinearBlock::_do_buildWorkSpace 0");
    4041    BysMcmc::Bsr::Gibbs::WorkSpace.ArimaReg wsar =
    4142    [[
    … …  
    4445      NameBlock _.arm.blk = arm.blk
    4546    ]];
     47  //WriteLn("TRACE"+_MID+" MainLinearBlock::_do_buildWorkSpace 1");
    4648    NameBlock _.ws := wsar;
     49  //WriteLn("TRACE"+_MID+" MainLinearBlock::_do_buildWorkSpace 2");
    4750    True
     51  };
     52
     53  ////////////////////////////////////////////////////////////////////////////
     54  Real _buildWorkSpace(VMatrix Y,    //Non filtered output
     55                       VMatrix X,    //Non filtered input
     56                       BysMcmc::BlockSampler arm.blk) //ArimaBlock filter
     57  ////////////////////////////////////////////////////////////////////////////
     58  {
     59  //WriteLn("TRACE"+_MID+" MainLinearBlock::_buildWorkSpace 0");
     60    If(!_.N1,
     61    {
     62    //WriteLn("TRACE"+_MID+" MainLinearBlock::_buildWorkSpace 1");
     63      _do_buildWorkSpace(Y,X,arm.blk)
     64    },
     65    {
     66    //WriteLn("TRACE"+_MID+" MainLinearBlock::_buildWorkSpace 2");
     67      VMatrix b1 = _.partialSampler::_.fixedValues;
     68      VMatrix X1 = SubCol(X, _.i1);
     69      VMatrix X2 = SubCol(X, _.i2);
     70      VMatrix Y2 = Y - X1*b1;
     71      _do_buildWorkSpace(Y2,X2,arm.blk)
     72    })
    4873  };
    4974  Code _.buildWorkSpace = _buildWorkSpace;

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

    ra2d7c82 r0c19991  
    2121  If(!Card(outputMissingBlock),
    2222  BysMcmc::Bsr::Gibbs::StdLinearBlock.Default(name, numBlock, firstCol, config),
    23   BysMcmc::Bsr::Gibbs::StdLinearFullBlock blk = [[
     23  BysMcmc::Bsr::Gibbs::StdLinearBlock blk = [[
    2424  //Given initial values of block variables
    2525  VMatrix _.store = Mat2VMat(SetMat(Extract(outputMissingBlock, 6)));

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

    ra2d7c82 r0c19991  
    103103  VMatrix _.beta = Rand(0,0,0,0);
    104104
     105  ////////////////////////////////////////////////////////////////////////////
    105106  Real putMID(Text mid)
     107  ////////////////////////////////////////////////////////////////////////////
    106108  {
    107109    Text _MID := mid;
    … …  
    120122  Real matchConstrains(VMatrix beta0, VMatrix A, VMatrix a)
    121123  ////////////////////////////////////////////////////////////////////////////
    122   {
    123    
     124  {   
    124125  //WriteLn("TRACE "+_MID+" matchConstrains 1 beta0="<<beta0);
    125126    Real ok = TRUE;
    … …  
    469470//////////////////////////////////////////////////////////////////////////////
    470471{
    471   BysMcmc::Bsr::Gibbs::StdLinearSubBlock subBlock =
     472  BysMcmc::Bsr::Gibbs::StdLinearSubBlock freeBlock =
    472473  [[ Text _MID = mid ]];
    473   subBlock
     474  freeBlock
    474475};
    475476
    … …  
    483484//////////////////////////////////////////////////////////////////////////////
    484485{
    485   Real initialize.subBlocks(Real unused);
    486   Real matchConstrains(VMatrix beta0, VMatrix A, VMatrix a);
    487 
    488   ////////////////////////////////////////////////////////////////////////////
    489   Real initialize(Real unused)
    490   ////////////////////////////////////////////////////////////////////////////
    491   {
    492     initialize.subBlocks(unused)
    493   };
    494486
    495487  ////////////////////////////////////////////////////////////////////////////
    … …  
    506498  {
    507499    _.ws
     500  };
     501
     502  BysMcmc::Bsr::Gibbs::StdLinearSubBlock _.fullBlock =
     503    BysMcmc::Bsr::Gibbs::StdLinearSubBlock.Build(_MID+"[_.fullBlock] ");
     504  BysMcmc::Bsr::Gibbs::StdLinearSubBlock _.freeBlock =
     505    BysMcmc::Bsr::Gibbs::StdLinearSubBlock.Build(_MID+"[_.freeBlock] ");
     506
     507  Set _.i1 = Copy(Empty); //_.partialSampler::getNonFreeIndex(0);
     508  Set _.i2 = Copy(Empty); //_.partialSampler::getFreeIndex(0);
     509  Real _.N1 = ?;
     510  Real _.N2 = ?;
     511
     512  ////////////////////////////////////////////////////////////////////////////
     513  Real initialize.subBlocks(Real unused)
     514  ////////////////////////////////////////////////////////////////////////////
     515  {
     516    Set _.i1 := _.partialSampler::getNonFreeIndex(0);
     517    Set _.i2 := _.partialSampler::getFreeIndex(0);
     518    Real _.N1 := Card(_.i1);
     519    Real _.N2 := Card(_.i2);
     520    Code _.freeBlock::standarize := standarize;
     521    Code _.freeBlock::getWorkSpace := getWorkSpace;
     522    Real _.freeBlock::usePreconditioning :=
     523      _.config::bsr.linBlk.usePreconditioning;
     524    Real _.freeBlock::cholesky.warningFreq :=
     525      _.config::bsr.cholesky.warningFreq;
     526    Real _.freeBlock::truncMNormal.gibbsNumIter :=
     527      _.config::bsr.truncMNormal.gibbsNumIter;
     528    Real _.freeBlock::drop_val :=
     529      _.config::bsr.linBlk.X.drop;
     530    Real _.freeBlock::_.eps :=
     531      _.config::bsr.cholesky.epsilon*_.N2*(_.N2+1);
     532    Real _.freeBlock::_.MH.useAfterIter :=
     533      Max(1,_.config::bsr.linBlk.MH.useAfterIter);
     534    Real _.freeBlock::_.MH.maxRejected :=
     535      _.config::bsr.linBlk.MH.maxRejected;
     536    Real _.freeBlock::_.MH.maxConsecutive :=
     537      _.config::bsr.linBlk.MH.maxConsecutive;
     538    1
     539  };
     540
     541  ////////////////////////////////////////////////////////////////////////////
     542  Real initialize(Real unused)
     543  ////////////////////////////////////////////////////////////////////////////
     544  {
     545    initialize.subBlocks(0)
     546  };
     547
     548  ////////////////////////////////////////////////////////////////////////////
     549  Real matchConstrains(VMatrix beta0, VMatrix A, VMatrix a)
     550  ////////////////////////////////////////////////////////////////////////////
     551  {
     552    If(!_.N1, 
     553    {
     554      Real ok = _.freeBlock::matchConstrains(beta0,A,a);
     555      Real putStoreV(_.freeBlock::_.beta0);
     556      ok
     557    },
     558    {
     559      Real ok1 = _.fullBlock::matchConstrains(beta0,A,a);
     560      Real ok2 = If(!ok1, False,
     561      {
     562        VMatrix b1 = _.partialSampler::_.fixedValues;
     563        VMatrix b2 = SubRow(_.fullBlock::_.beta0, _.i2);
     564        VMatrix A1 = SubCol(A, _.i1);
     565        VMatrix A2 = SubCol(A, _.i2);
     566        VMatrix a2 = a - A1*b1;
     567        Real _.freeBlock::matchConstrains(b2,A2,a2)
     568      });
     569      VMatrix b_ = _.partialSampler::mergeWithFreeSample(_.freeBlock::_.beta0);
     570      Real putStoreV(b_);
     571      ok2
     572    })
     573  };
     574
     575  ////////////////////////////////////////////////////////////////////////////
     576  VMatrix draw(Real numSim)
     577  ////////////////////////////////////////////////////////////////////////////
     578  {
     579  //WriteLn("TRACE"+_MID+" StdLinearBlock::draw 0");
     580    If(!_.N1, 
     581    {
     582    //WriteLn("TRACE"+_MID+" StdLinearBlock::draw 1.1");
     583    _.freeBlock::draw(numSim)
     584    },
     585    {
     586    //WriteLn("TRACE"+_MID+" StdLinearBlock::draw 2.1");
     587      VMatrix b1 = _.partialSampler::_.fixedValues;
     588    //WriteLn("TRACE"+_MID+" StdLinearBlock::draw 2.2");
     589      VMatrix A1 = SubCol(_.fullBlock::_.A, _.i1);
     590    //WriteLn("TRACE"+_MID+" StdLinearBlock::draw 2.3");
     591      VMatrix _.freeBlock::_.a := _.fullBlock::_.a - A1*b1;
     592    //WriteLn("TRACE"+_MID+" StdLinearBlock::draw 2.4");
     593      _.freeBlock::draw(numSim)
     594    })
    508595  }
    509596
    510597};
     598
    511599
    512600//////////////////////////////////////////////////////////////////////////////
    … …  
    529617  aux
    530618};
    531 
    532 //////////////////////////////////////////////////////////////////////////////
    533 Text _.autodoc.member.StdLinearFullBlock =
    534 "Builds a NameBlock that can draw a Gibbs sample of an standarized linear "
    535 "block of a general regression model with variant input and ouput regression "
    536 "matrices and optional but fixed constraining inequations.";
    537 //////////////////////////////////////////////////////////////////////////////
    538 Class StdLinearFullBlock: StdLinearBlock
    539 //////////////////////////////////////////////////////////////////////////////
    540 {
    541   BysMcmc::Bsr::Gibbs::StdLinearSubBlock _.subBlock =
    542     BysMcmc::Bsr::Gibbs::StdLinearSubBlock.Build(_MID+"[_.subBlock] ");
    543 
    544   ////////////////////////////////////////////////////////////////////////////
    545   //Default behaviour assumes standarized linear regression
    546   VMatrix standarize(VMatrix _Z)
    547   ////////////////////////////////////////////////////////////////////////////
    548   {
    549     _Z
    550   };
    551 
    552   ////////////////////////////////////////////////////////////////////////////
    553   NameBlock getWorkSpace(Real unused)
    554   ////////////////////////////////////////////////////////////////////////////
    555   {
    556     _.ws
    557   };
    558 
    559   ////////////////////////////////////////////////////////////////////////////
    560   Real initialize.subBlocks(Real unused)
    561   ////////////////////////////////////////////////////////////////////////////
    562   {
    563     Real _.subBlock::putMID(_MID+"[_.subBlock] ");
    564     Code _.subBlock::standarize := standarize;
    565     Code _.subBlock::getWorkSpace := getWorkSpace;
    566     Real _.subBlock::usePreconditioning :=
    567       _.config::bsr.linBlk.usePreconditioning;
    568     Real _.subBlock::cholesky.warningFreq :=
    569       _.config::bsr.cholesky.warningFreq;
    570     Real _.subBlock::truncMNormal.gibbsNumIter :=
    571       _.config::bsr.truncMNormal.gibbsNumIter;
    572     Real _.subBlock::drop_val :=
    573       _.config::bsr.linBlk.X.drop;
    574     Real _.subBlock::_.eps :=
    575       _.config::bsr.cholesky.epsilon*_.N*(_.N+1);
    576     Real _.subBlock::_.MH.useAfterIter :=
    577       Max(1,_.config::bsr.linBlk.MH.useAfterIter);
    578     Real _.subBlock::_.MH.maxRejected :=
    579       _.config::bsr.linBlk.MH.maxRejected;
    580     Real _.subBlock::_.MH.maxConsecutive :=
    581       _.config::bsr.linBlk.MH.maxConsecutive;
    582     1
    583   };
    584 
    585   ////////////////////////////////////////////////////////////////////////////
    586   Real initialize(Real unused)
    587   ////////////////////////////////////////////////////////////////////////////
    588   {
    589     initialize.subBlocks(unused)
    590   };
    591 
    592   ////////////////////////////////////////////////////////////////////////////
    593   Real matchConstrains(VMatrix beta0, VMatrix A, VMatrix a)
    594   ////////////////////////////////////////////////////////////////////////////
    595   {
    596     Real _.subBlock::matchConstrains(beta0,A,a);
    597     Real putStoreV(_.subBlock::_.beta0);
    598     1
    599   };
    600 
    601   ////////////////////////////////////////////////////////////////////////////
    602   VMatrix draw(Real numSim)
    603   ////////////////////////////////////////////////////////////////////////////
    604   {
    605     _.subBlock::draw(numSim)
    606   }
    607 
    608 };
    609 /*
    610 //////////////////////////////////////////////////////////////////////////////
    611 Text _.autodoc.member.StdLinearPartialBlock =
    612 "Builds a NameBlock that can draw a Gibbs sample of an standarized linear "
    613 "block of a general regression model that is partially fixed or reloaded.";
    614 //////////////////////////////////////////////////////////////////////////////
    615    Class StdLinearPartialBlock: BlockSampler
    616 //////////////////////////////////////////////////////////////////////////////
    617 {
    618   BysMcmc::StdLinearSubBlock _.subBlock = StdLinearSubBlock.Build(0);
    619   B
    620 };
    621 */
Note: See TracChangeset for help on using the changeset viewer.