Changeset c9971f2 in tol


Ignore:
Timestamp:
Oct 23, 2009, 7:45:26 AM (17 years ago)
Author:
Víctor de Buen <vdebuen@…>
Branches:
1941-docker-image-to-test, devel, master
Children:
e9c0ec3
Parents:
16b4825
Message:

Preparing for joint input and output missing simulation
Related to #759

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

Location:
tol/stdlib/math/stat/models/bayesian/bysMcmc/bsr/gibbs
Files:
2 deleted
5 edited

Legend:

Unmodified
Added
Removed

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

    r16b4825 rc9971f2  
    2121  BysMcmc::Bsr::Gibbs::@BlockSigma  _.sig.blk;
    2222  BysMcmc::Bsr::Gibbs::@BlockArima  _.arm.blk;
    23   BysMcmc::Bsr::Gibbs::@BlockStdLin _.inpMis.blk;
    24   BysMcmc::Bsr::Gibbs::@BlockStdLin _.outMis.blk;
     23  BysMcmc::Bsr::Gibbs::@BlockStdLin _.mis.blk;
    2524  //Non linear blocks
    2625  Set _.nonLinBlk;
    … …  
    157156      });
    158157      Set nonLinFilter = BinGroup("<<", aux.nlb);
    159       Real nextCol = _.outMis.blk::_.lastCol+1;
     158      Real nextCol = _.mis.blk::_.lastCol+1;
    160159      @MasterSamplerBasic T = _this;
    161160      For(1,Card(nonLinFilter), NameBlock(Real numBlk)
    … …  
    186185      _.sig.blk,
    187186      _.arm.blk,
    188       _.inpMis.blk,
    189       _.outMis.blk
     187      _.mis.blk
    190188    ]]
    191189    << _.nonLinBlk;
    … …  
    213211
    214212  //WriteLn("TRACE"+_MID+" initialize 1");
    215     VMatrix _.YplusV := _.outMis.blk::filter(_.Y);
     213    VMatrix _.YplusV := _.mis.blk::filter.output(_.Y);
    216214  //WriteLn("TRACE"+_MID+" initialize 1.1 _.YplusV("<<VRows(_.YplusV)+"x"<<VColumns(_.YplusV)+")");
    217     VMatrix _.XplusU := _.inpMis.blk::filter(_.modelDef->X);
     215    VMatrix _.XplusU := _.mis.blk::filter.input(_.modelDef->X);
    218216  //WriteLn("TRACE"+_MID+" initialize 1.2 _.XplusU("<<VRows(_.XplusU)+"x"<<VColumns(_.XplusU)+")");
    219217    Matrix givenBeta=
    … …  
    267265  //WriteLn("TRACE"+_MID+" resume 0");
    268266    Real initialize.blocks(unused);
    269     VMatrix _.YplusV := _.outMis.blk::filter(_.Y);
    270     VMatrix _.XplusU := _.inpMis.blk::filter(_.modelDef->X);
     267    VMatrix _.YplusV := _.mis.blk::filter.output(_.Y);
     268    VMatrix _.XplusU := _.mis.blk::filter.input(_.modelDef->X);
    271269    VMatrix _.beta   := _.mln.blk::_.store;
    272270    Real ok = _.mln.blk::matchConstrains(_.beta, _.modelDef->A, _.modelDef->a);
    … …  
    305303  //WriteLn("TRACE"+_MID+" @MasterSamplerBasic::draw 1");
    306304    Real DoNonLinOutputFilter(0);
    307     VMatrix om = If(Not(_.outMis.blk::_.N),
     305
     306    VMatrix mis = If(Not(_.mis.blk::_.N),
    308307    {
    309308      VMatrix _.YplusV := _.Y;
    310       Rand(0,1,0,0)
    311     },
    312     {
    313       Real If(_.outMis.blk::_.partialSampler::_.freeParamNum,
    314       {
    315         Real _.outMis.blk::_.buildWorkSpace(
    316           _.arm.blk,_.sig.blk::_.sigmaInv,_.Y-_.XplusU*_.beta)
    317       });
    318       VMatrix aux = _.outMis.blk::get.draw(numSim, False);
    319       VMatrix _.YplusV := _.outMis.blk::filter(_.Y);
    320       aux
    321     });
    322   //WriteLn("TRACE"+_MID+" @MasterSamplerBasic::draw 2 om="<<om);
    323     VMatrix im = If(Not(_.inpMis.blk::_.N),
    324     {
    325309      VMatrix _.XplusU := _.modelDef->X;
    326310      Rand(0,1,0,0)
    327311    },
    328312    {
    329       Real If(_.inpMis.blk::_.partialSampler::_.freeParamNum,
    330       {
    331         Real _.inpMis.blk::_.buildWorkSpace(
    332           _.arm.blk,_.sig.blk::_.sigmaInv,_.YplusV-_.modelDef->X*_.beta,_.beta)
    333       });
    334       VMatrix aux = _.inpMis.blk::get.draw(numSim, False);
    335       VMatrix _.XplusU := _.inpMis.blk::filter(_.modelDef->X);
     313      Real If(_.mis.blk::_.partialSampler::_.freeParamNum,
     314      {
     315        Real _.mis.blk::_.buildWorkSpace(
     316          _.arm.blk,_.sig.blk::_.sigmaInv,_.modelDef->Y-_.modelDef->X*_.beta,_.beta)
     317      });
     318      VMatrix aux = _.mis.blk::get.draw(numSim, False);
     319      VMatrix _.YplusV := _.mis.blk::filter.output(_.Y);
     320      VMatrix _.XplusU := _.mis.blk::filter.input(_.modelDef->X);
    336321      aux
    337322    });
     323
    338324  //WriteLn("TRACE"+_MID+" @MasterSamplerBasic::draw 3 im="<<im);
    339325    VMatrix _.filter := If(Not(_.mln.blk::_.N),
    … …  
    390376    });
    391377  //WriteLn("TRACE"+_MID+" @MasterSamplerBasic::draw 6 _.beta="<<_.beta);
    392     VMatrix param = _.beta << sig << arm << im << om << nlb;
     378    VMatrix param = _.beta << sig << arm << mis << nlb;
    393379  //WriteLn("TRACE"+_MID+" @MasterSamplerBasic::draw 7 param="<<VRows(param)+"x"<<VColumns(param)+")");
    394380    Real putStoreV(param);
    … …  
    428414    Real DoNonLinOutputFilter(0);
    429415  //WriteLn("TRACE"+_MID+" eval 2");
    430     VMatrix _.YplusV := _.outMis.blk::filter(_.Y);
    431     VMatrix fullOutput.mis = _.outMis.blk::filter(_.Y*0);
     416    VMatrix _.YplusV := _.mis.blk::filter.output(_.Y);
     417    VMatrix fullOutput.mis = _.mis.blk::filter.output(_.Y*0);
    432418  //WriteLn("TRACE"+_MID+" eval 3");
    433     VMatrix _.XplusU := _.inpMis.blk::filter(_.modelDef->X);
     419    VMatrix _.XplusU := _.mis.blk::filter.input(_.modelDef->X);
    434420  //WriteLn("TRACE"+_MID+" eval 4");
    435421    VMatrix _.beta := _.mln.blk::_.store;
    … …  
    516502          Serie serResiduals = MatSerSet(Tra(residuals),ti->Dating, difFirst)[1];
    517503          Serie serForecast  = MatSerSet(Tra(forecast), ti->Dating, difFirst)[1];
    518           [[sigma, stdErr, avr, serOutput, serOutput.mis, serOutput.full,
     504          [[sigma, stdErr, avr, serOutput, serOutput.mis, serOutput.extFil, serOutput.full,
    519505            serFilter, serFilter.mis, serFilter.full, serNoise, serForecast,
    520506            serDifNoise, serResiduals, arima]]
    … …  
    694680  ////////////////////////////////////////////////////////////////////////////
    695681  {
    696     VMatrix X_U = _.inpMis.blk::filter(_.modelDef->X);
     682    VMatrix X_U = _.mis.blk::filter.input(_.modelDef->X);
    697683    VMatrix X = SubCol(X_U,beta.index);
    698684    VMat2Mat(subRegEquBlockV(resBlock, X))
    … …  
    707693  ////////////////////////////////////////////////////////////////////////////
    708694  {
    709     VMatrix Y = _.outMis.blk::filter(_.modelDef->Y);
     695    VMatrix Y = _.mis.blk::filter.output(_.modelDef->Y);
    710696    VMat2Mat(subRegEquBlockV(resBlock, Y))
    711697  };
    … …  
    733719  ////////////////////////////////////////////////////////////////////////////
    734720  {
    735     VMatrix X_U = _.inpMis.blk::filter(_.modelDef->X);
     721    VMatrix X_U = _.mis.blk::filter.input(_.modelDef->X);
    736722    VMat2Mat(subDifRegEquBlockV(resBlock, X_U))
    737723  };
    … …  
    745731  ////////////////////////////////////////////////////////////////////////////
    746732  {
    747     VMatrix X_U = _.inpMis.blk::filter(_.modelDef->X);
     733    VMatrix X_U = _.mis.blk::filter.input(_.modelDef->X);
    748734    VMatrix X = SubCol(X_U,beta.index);
    749735    VMat2Mat(subDifRegEquBlockV(resBlock, X))
    … …  
    758744  ////////////////////////////////////////////////////////////////////////////
    759745  {
    760     VMatrix Y_V = _.outMis.blk::filter(_.modelDef->Y);
     746    VMatrix Y_V = _.mis.blk::filter.output(_.modelDef->Y);
    761747    VMat2Mat(subDifRegEquBlockV(resBlock, Y_V))
    762748  };
    … …  
    898884      );
    899885
    900       BysMcmc::Bsr::Gibbs::@BlockStdLin _.inpMis.blk =
    901       BysMcmc::Bsr::Gibbs::@BlockInputMissing::New
     886      BysMcmc::Bsr::Gibbs::@BlockStdLin _.mis.blk =
     887      BysMcmc::Bsr::Gibbs::@BlockMissing::New
    902888      (
    903         "@BSR.InputMissingBlock."+name,
     889        "@BSR.MissingBlock."+name,
     890        _.modelDef->OutputMissingBlock,
    904891        _.modelDef->InputMissingBlock,
     892        VRows(_.modelDef->Y),
    905893        mcmcPartialHandler,
    906894        4,
    907895        _.arm.blk::_.lastCol+1,
    908         config
    909       );
    910 
    911       BysMcmc::Bsr::Gibbs::@BlockStdLin _.outMis.blk =
    912       BysMcmc::Bsr::Gibbs::@BlockOutputMissing::New
    913       (
    914         "@BSR.OutputMissingBlock."+name,
    915         _.modelDef->OutputMissingBlock,
    916         mcmcPartialHandler,
    917         5,
    918         _.inpMis.blk::_.lastCol+1,
    919896        config
    920897      );

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

    r16b4825 rc9971f2  
    1212  "Gibbs related methods for BSR and extensions.";
    1313
     14  //Sigma block related methods
    1415  #Embed "_sigmaBlock.tol";
    1516  //ARIMA block related methods
    … …  
    1920  //Linear blocks related methods
    2021  #Embed "_mainLinearBlock.tol";
    21   //Sigma block related methods
    22   //Missing input block related methods
    23   #Embed "_inputMissingBlock.tol";
    24   //Missing output block related methods
    25   #Embed "_outputMissingBlock.tol";
    26   //Generic I/O non linear filter block
     22  //Missing block related methods
     23  #Embed "_inputMissing.tol";
     24  #Embed "_outputMissing.tol";
     25  #Embed "_missingBlock.tol";
     26  //Generic non linear filter block
    2727  #Embed "_nonLinBlock.tol";
    2828  //Basic BSR master (Lineal+Sigma+ARIMA+Missing I/O)

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

    r16b4825 rc9971f2  
    179179  ////////////////////////////////////////////////////////////////////////////
    180180  {
    181     VMatrix _.H := arm.blk::filter(-si*_.HEye);
     181    VMatrix _.H := si*arm.blk::filter(_.HEye);
    182182    True
    183183  };
    … …  
    197197    VMatrix _.H := If(!_.n2, Zeros(_.M, _.N2),
    198198    {
    199       arm.blk::filter(si*_.HEye)
     199      si*arm.blk::filter(_.HEye)
    200200    });
    201201    True

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

    r16b4825 rc9971f2  
    233233    If(_.N<=0,0,{
    234234    Real putStore(values);
    235     VMatrix sigmaInvVar = BinGroup("+",
    236     For(1,Card(_.sigmaBlockVar), VMatrix (Real k)
     235    VMatrix sigmaInvVar = SetSum(For(1,Card(_.sigmaBlockVar), VMatrix (Real k)
    237236    {
    238237      Real si = 1/Sqrt(VMatDat(_.store,k,1));
    … …  
    292291    //regression standarizing
    293292    VMatrix _.sigmaInvCte = If(!Card(_.sigmaBlockCte), Zeros(_.m),
    294       BinGroup("+", EvalSet(_.sigmaBlockCte, VMatrix (@BSR.NoiseDistrib d)
     293      SetSum(EvalSet(_.sigmaBlockCte, VMatrix (@BSR.NoiseDistrib d)
    295294      { d->Li })));
    296295    //Sigma block current store

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

    r16b4825 rc9971f2  
    592592  Real _.N2 = ?;
    593593
     594  //Indexes of inequations related to free variables
     595  Set _.r2 = Copy(Empty);
     596
    594597  //--------------------------------------------------------------------------
    595598  //Methods
    … …  
    666669      Real ok2 = If(!ok1, False,
    667670      {
    668         If(!_.N2, True,
     671        Case(!_.N2, True,
     672        !VRows(A),
     673        {
     674          VMatrix b2 = SubRow(_.fullBlock::_.beta0, _.i2);
     675          _.freeBlock::matchConstrains(b2,Rand(0,_.N2,0,0),Rand(0,1,0,0))
     676        },
     677        1==1, 
    669678        {
    670679        //WriteLn("TRACE"+_MID+" @BlockStdLin::matchConstrains 2.4");
    … …  
    673682          VMatrix A1 = SubCol(A, _.i1);
    674683        //WriteLn("TRACE"+_MID+" @BlockStdLin::matchConstrains 2.6");
    675           VMatrix A2 = SubCol(A, _.i2);
     684          VMatrix A2_ = SubCol(A, _.i2);
     685          VMatrix a2_ = a - A1*b1;
    676686        //WriteLn("TRACE"+_MID+" @BlockStdLin::matchConstrains 2.7");
    677           VMatrix a2 = a - A1*b1;
     687          Matrix aux = VMat2Mat(Abs(A2_)*Rand(_.N2,1,1,1));
     688          Matrix nz  = NE(aux,Rand(_.fullBlock::_.R,1,0,0));
     689          Matrix idx = SetCol(Range(1,_.fullBlock::_.R,1));
     690          Set nz_idx = MatSet(Tra(nz $* idx))[1];
     691          Set _.r2 = Select(nz_idx,Real(Real j) { j>0 });
     692          VMatrix A2 = SubRow(A2_,_.r2);
     693          VMatrix a2 = SubRow(a2_,_.r2);
    678694        //WriteLn("TRACE"+_MID+" @BlockStdLin::matchConstrains 2.8");
    679           Real _.freeBlock::matchConstrains(b2,A2,a2)
     695          _.freeBlock::matchConstrains(b2,A2,a2)
    680696       })
    681697      });
    … …  
    705721      //WriteLn("TRACE"+_MID+" @BlockStdLin::draw 2.1");
    706722        VMatrix b1 = SubRow(_.partialSampler::_.values, _.i1);
    707       //WriteLn("TRACE"+_MID+" @BlockStdLin::draw 2.2");
    708         VMatrix A1 = SubCol(_.fullBlock::_.A, _.i1);
    709       //WriteLn("TRACE"+_MID+" @BlockStdLin::draw 2.3");
    710         VMatrix _.freeBlock::_.a := _.fullBlock::_.a - A1*b1;
     723        Real If(_.freeBlock::_.R,
     724        {
     725        //WriteLn("TRACE"+_MID+" @BlockStdLin::draw 2.2");
     726          VMatrix A1 = SubCol(_.fullBlock::_.A, _.i1);
     727        //WriteLn("TRACE"+_MID+" @BlockStdLin::draw 2.3");
     728          VMatrix a2_ = _.fullBlock::_.a - A1*b1;
     729          VMatrix _.freeBlock::_.a := SubRow(a2_,_.r2);
     730          True
     731        });
    711732      //WriteLn("TRACE"+_MID+" @BlockStdLin::draw 2.4");
    712733        True
Note: See TracChangeset for help on using the changeset viewer.