Changeset b3d4b55 in tol


Ignore:
Timestamp:
Oct 14, 2009, 2:45:47 PM (17 years ago)
Author:
Víctor de Buen <vdebuen@…>
Branches:
1941-docker-image-to-test, devel, master
Children:
dd61b24
Parents:
1904f20
Message:

Unifying common issues of input and output missing handling

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

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

Legend:

Unmodified
Added
Removed

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

    r1904f20 rb3d4b55  
    1515  //Set of @Bsr.Missing.Info without repeatitions
    1616  Set _.missingInfo;
     17  Real _.prior.M;
     18  //Average vector of missing priors
     19  Matrix _.prior.nu;
     20  //Inverse of standard deviation vector of missing priors
     21  Matrix _.prior.si;
     22  //Standarized output of missing priors
     23  VMatrix _.prior.Y;
     24  //Standarized input of missing priors
     25  VMatrix _.prior.X;
    1726  //Classification of missing by name
    1827  Set _.inputMissingDistinct;
    … …  
    2837  Set _.prior.normal;
    2938  //Number of input missing prior equations
    30   Real _.prior.M;
    31   //Average vector of input missing priors
    32   Matrix _.prior.nu;
    33   //Inverse of standard deviation vector of input missing priors
    34   Matrix _.prior.si;
    35   //Standarized output of input missing priors
    36   VMatrix _.prior.Y;
    37   //Standarized input of input missing priors
    38   VMatrix _.prior.X;
    39 
    4039
    4140  //--------------------------------------------------------------------------
    4241  // Auxiliar members
    4342  //--------------------------------------------------------------------------
     43  //Work space
    4444  @WorkSpace.LinReg _.ws = BysMcmc::Bsr::Gibbs::@WorkSpace.LinReg::Null(0);
    4545  //Standarized input of free missing priors
    … …  
    148148  //WriteLn("TRACE "+_MID+" 7");
    149149    VMatrix store = _.store;
    150     //Selecting input missing variables with left truncated prior information
     150    //Selecting missing variables with left truncated prior information
    151151    Set leftConstrained =
    152152      Select(_.missingInfo, Real(@BSR.MissingBlock b)
    … …  
    161161      });
    162162  //WriteLn("TRACE "+_MID+" initialize 2");
    163     //Number of input missing with left constrained prior
     163    //Number of missing with left constrained prior
    164164    Real left.R = Card(leftConstrained);
    165     //Number of input missing with left constrained prior
     165    //Number of missing with left constrained prior
    166166    Real right.R = Card(rightConstrained);
    167     //Number of input missing with constrained prior
     167    //Number of missing with constrained prior
    168168    Real R = left.R + right.R;
    169     //Border of input missing constrains
     169    //Border of missing constrains
    170170    VMatrix prior.a = If(!R, Rand(0,0,0,0),
    171171    {
    … …  
    179179    });
    180180  //WriteLn("TRACE "+_MID+" initialize 3");
    181     //Coefficients of input missing constrains
     181    //Coefficients of missing constrains
    182182    VMatrix prior.A = If(!R, Rand(0,0,0,0),
    183183    {
    … …  
    204204  ////////////////////////////////////////////////////////////////////////////
    205205  Text _.autodoc.member.filter =
    206   "Filter regression input with current missing block values. You can set a "
     206  "Filter regression with current missing block values. You can set a "
    207207  "previously simulated vector of paramters calling before to method "
    208208  "setMcmcRow";
    … …  
    210210  ////////////////////////////////////////////////////////////////////////////
    211211  {
    212   //WriteLn("TRACE "+_MID+" filter 1");
     212  //WriteLn("TRACE"+_MID+" filter Z:"<<Z);
    213213    If(Not(enabled), Z,
    214214    {
    … …  
    221221  };
    222222
    223   //////////////////////////////////////////////////////////////////////////////
    224   Static Text _.autodoc.member.New =
     223  ////////////////////////////////////////////////////////////////////////////
     224  Text _.autodoc.member.New =
    225225  "Builds a NameBlock that can draw a Gibbs sample of a Input Missing block "
    226226  "of a Bayesian Sparse Regression model";
    … …  
    238238    BysMcmc::Bsr::Gibbs::@BlockStdLinEmpty::New(name, numBlock, firstCol, config),
    239239    BysMcmc::Bsr::Gibbs::@BlockInputMissing blk = [[
    240     Code _.buildWorkSpace = _buildWorkSpace;
    241240    Set _.inputMissingDistinct =
    242241    {
    … …  
    265264    Set _.storeIndex = EvalSet(inputMissingBlock,
    266265     Real(@BSR.MissingBlock mis){ mis->Index }),
    267     //Given initial values of block variables
    268     VMatrix _.store = Mat2VMat(SetMat(
    269        EvalSet(_.inputMissingDistinct, Set(Set cl){ [[(cl[1])->Nu]] })));
    270     //Defining the block in the Markov Chain
    271     Real _defined = define(
    272       name,
    273       EvalSet(_.inputMissingDistinct, Text(Set cl){ (cl[1])->Name }),
    274       mcmcPartialHandler,
    275       numBlock,
    276       firstCol,
    277       config);
    278266    //Input missing pairs (equation,variable)
    279267    Matrix _.tj = {
    … …  
    289277    //Input missing index j
    290278    Set _.j = Traspose(Extract(inputMissingBlock, 4))[1];
    291 
    292       Set _.missingInfo = EvalSet(_.inputMissingDistinct, Set(Set cl){ cl[1] });
     279    Set _.missingInfo = EvalSet(_.inputMissingDistinct, Set(Set cl){ cl[1] });
     280
     281    //Given initial values of block variables
     282    VMatrix _.store = Mat2VMat(SetMat(Extract(_.missingInfo, 6)));
     283    Code _.buildWorkSpace = _buildWorkSpace;
     284    //Defining the block in the Markov Chain
     285    Real _defined = define(
     286      name,
     287      EvalSet(_.missingInfo, Text(@BSR.MissingBlock b){ b->Name }),
     288      mcmcPartialHandler,
     289      numBlock,
     290      firstCol,
     291      config);
    293292    //Selecting input missing variables with prior information
    294     Set _.prior.normal =
    295     {
    296       Select(_.missingInfo, Real(@BSR.MissingBlock b)
    297       {
    298         Or(b->Prior == "Normal",b->Prior == "TruncatedNormal")
    299       })
    300     };
    301     //Number of input missing prior equations
     293    Set _.prior.normal = Select(_.missingInfo, Real(@BSR.MissingBlock b)
     294    {
     295      Or(b->Prior == "Normal",b->Prior == "TruncatedNormal")
     296    });
     297    //Number of missing prior equations
    302298    Real _.prior.M = Card(_.prior.normal);
    303     //Average vector of input missing priors
    304     Matrix _.prior.nu = If(Not(_.prior.M),Rand(0,1,0,0),
     299    //Average vector of missing priors
     300    Matrix _.prior.nu = If(!_.prior.M,Rand(0,1,0,0),
    305301    {
    306302    //WriteLn("TRACE "+_MID+" 3");
    307303      SetMat(Extract(_.prior.normal, 6))
    308304    });
    309     //Inverse of standard deviation vector of input missing priors
    310     Matrix _.prior.si = If(Not(_.prior.M),Rand(0,1,0,0),
     305    //Inverse of standard deviation vector of missing priors
     306    Matrix _.prior.si = If(!_.prior.M,Rand(0,1,0,0),
    311307    {
    312308    //WriteLn("TRACE "+_MID+" 4");
    … …  
    314310      RPow(s2,-1/2)
    315311    });
    316     //Standarized output of input missing priors
    317     VMatrix _.prior.Y = If(Not(_.prior.M),Rand(0,1,0,0),
     312    //Standarized output of missing priors
     313    VMatrix _.prior.Y = If(!_.prior.M,Rand(0,1,0,0),
    318314      Mat2VMat(_.prior.nu $* _.prior.si));
    319     //Standarized input of input missing priors
    320     VMatrix _.prior.X = If(Not(_.prior.M),Rand(0,_.N,0,0),
     315    //Standarized input of missing priors
     316    VMatrix _.prior.X = If(!_.prior.M,Rand(0,_.N,0,0),
    321317    {
    322318    //WriteLn("TRACE "+_MID+" _.prior.X 1");
    323     //WriteLn("TRACE "+_MID+" 5");
    324       Matrix aux1= SetMat(Traspose([[Set Range(1,_.prior.M,1)]]));
    325     //WriteLn("TRACE "+_MID+" 6");
    326       Matrix aux2= SetMat(Extract(_.prior.normal, 2));
    327       Matrix triplet =  aux1 | aux2 | _.prior.si;
     319      Matrix triplet = SetMat(Traspose([[Set Range(1,_.prior.M,1)]])) |
     320                       SetMat(Extract(_.prior.normal, 2)) |
     321                       _.prior.si;
    328322    //WriteLn("TRACE "+_MID+" _.prior.X 2 triplet=\n"<<triplet);
    329323      VMatrix tr = Triplet(triplet,_.prior.M,_.N);

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

    r1904f20 rb3d4b55  
    1313  // Mandatory members
    1414  //--------------------------------------------------------------------------
    15   //Set of @Bsr.Missing.Info
     15  //Set of @Bsr.Missing.Info without repeatitions
    1616  Set _.missingInfo;
     17  Real _.prior.M;
     18  //Average vector of missing priors
     19  Matrix _.prior.nu;
     20  //Inverse of standard deviation vector of missing priors
     21  Matrix _.prior.si;
     22  //Standarized output of missing priors
     23  VMatrix _.prior.Y;
     24  //Standarized input of missing priors
     25  VMatrix _.prior.X;
    1726  //Output missing pairs i,k
    1827  Set _.t;
    … …  
    2029  Matrix _.I;
    2130  //Number of output missing prior equations
    22   Real _.prior.M;
    23   //Average vector of output missing priors
    24   Matrix _.prior.nu;
    25   //Inverse of standard deviation vector of output missing priors
    26   Matrix _.prior.si;
    27   //Standarized output of missing priors
    28   VMatrix _.prior.Y;
    29   //Standarized input of missing priors
    30   VMatrix _.prior.X;
    31 
    32   //--------------------------------------------------------------------------
    33   //Auxiliar members
    34   //--------------------------------------------------------------------------
    35   //Auxiliar variable for filter creating
    36   VMatrix _.FEye = Rand(0,0,0,0);
     31
     32  //--------------------------------------------------------------------------
     33  // Auxiliar members
     34  //--------------------------------------------------------------------------
    3735  //Work space
    3836  @WorkSpace.LinReg _.ws = BysMcmc::Bsr::Gibbs::@WorkSpace.LinReg::Null(0);
    … …  
    4038  VMatrix _.prior.X.free = Rand(0,0,0,0);
    4139
     40  //Auxiliar variable for filter creating
     41  VMatrix _.FEye = Rand(0,0,0,0);
    4242  //--------------------------------------------------------------------------
    4343  //Methods
    … …  
    125125    Real initialize.subBlocks(unused);
    126126  //WriteLn("TRACE "+_MID+" initialize 1");
    127     VMatrix u = _.store;
    128     //Selecting output missing variables with left truncated prior information
     127  //WriteLn("TRACE "+_MID+" 7");
     128    VMatrix store = _.store;
     129    //Selecting missing variables with left truncated prior information
    129130    Set leftConstrained =
    130131      Select(_.missingInfo, Real(@BSR.MissingBlock b)
    … …  
    139140      });
    140141  //WriteLn("TRACE "+_MID+" initialize 2");
    141     //Number of output missing with left constrained prior
     142    //Number of missing with left constrained prior
    142143    Real left.R = Card(leftConstrained);
    143     //Number of output missing with left constrained prior
     144    //Number of missing with left constrained prior
    144145    Real right.R = Card(rightConstrained);
    145     //Number of output missing with constrained prior
     146    //Number of missing with constrained prior
    146147    Real R = left.R + right.R;
    147     //Border of output missing constrains
     148    //Border of missing constrains
    148149    VMatrix prior.a = If(!R, Rand(0,0,0,0),
    149150    {
    150       If(!left.R,  Rand(0,1,0,0), (-Mat2VMat(SetMat(Extract(leftConstrained,  8)))))  <<
    151       If(!right.R, Rand(0,1,0,0), (+Mat2VMat(SetMat(Extract(rightConstrained, 9))))) 
     151    //WriteLn("TRACE "+_MID+" 8");
     152      VMatrix aux1=If(!left.R,  Rand(0,1,0,0),
     153        (-Mat2VMat(SetMat(Extract(leftConstrained,  8)))));
     154    //WriteLn("TRACE "+_MID+" 9");
     155      VMatrix aux2=If(!right.R, Rand(0,1,0,0),
     156        (+Mat2VMat(SetMat(Extract(rightConstrained, 9)))));
     157      aux1<<aux2
    152158    });
    153159  //WriteLn("TRACE "+_MID+" initialize 3");
    154     //Coefficients of output missing constrains
    155     VMatrix prior.A =
     160    //Coefficients of missing constrains
     161    VMatrix prior.A = If(!R, Rand(0,0,0,0),
    156162    {
    157163      Set left.range  = Range(1,        left.R, 1);
    158164      Set right.range = Range(1+left.R, R,      1);
    159       Matrix left  = If(!left.R, Rand(0,3,0,0),
    160                      SetMat(Traspose([[left.range]])) |
    161                      SetMat(Extract(leftConstrained, 2)) |
    162                      Rand(left.R,1,-1,-1));
     165    //WriteLn("TRACE "+_MID+" 10");
     166      Matrix left  = If(!left.R, Rand(0,3,0,0),
     167                                 SetMat(Traspose([[left.range]])) |
     168                                 SetMat(Extract(leftConstrained, 2)) |
     169                                 Rand(left.R,1,-1,-1));
     170    //WriteLn("TRACE "+_MID+" 11");
    163171      Matrix right = If(!right.R, Rand(0,3,0,0),
    164                      SetMat(Traspose([[right.range]])) |
    165                      SetMat(Extract(rightConstrained, 2)) |
    166                      Rand(right.R,1,+1,+1));
     172                                  SetMat(Traspose([[right.range]])) |
     173                                  SetMat(Extract(rightConstrained, 2)) |
     174                                  Rand(right.R,1,+1,+1));
    167175      Convert(Triplet(left<<right,R,_.N),"Cholmod.R.Sparse")
    168     };
     176    });
    169177  //WriteLn("TRACE "+_MID+" initialize 4");
    170     Real ok = matchConstrains(u,prior.A,prior.a);
     178    Real ok = matchConstrains(store,prior.A,prior.a);
    171179  //WriteLn("TRACE "+_MID+" initialize 5");
    172180    ok
    173181  };
    174  
     182
    175183  ////////////////////////////////////////////////////////////////////////////
    176184  Text _.autodoc.member.filter =
    177   "Filter regression output with current missing block values. You can set a "
     185  "Filter regression with current missing block values. You can set a "
    178186  "previously simulated vector of paramters calling before to method "
    179187  "setMcmcRow";
    … …  
    193201  "Builds a NameBlock that can draw a Gibbs sample of a Output Missing block "
    194202  "of a Bayesian Sparse Regression model";
    195   ////////////////////////////////////////////////////////////////////////////
     203  //////////////////////////////////////////////////////////////////////////////
    196204  Static @BlockStdLin New(
    197205    Text name,
    … …  
    206214    BysMcmc::Bsr::Gibbs::@BlockStdLinEmpty::New(name, numBlock, firstCol, config),
    207215    BysMcmc::Bsr::Gibbs::@BlockOutputMissing blk = [[
     216    //Output missing pairs i,k
     217    Set _.t = Traspose(Extract(outputMissingBlock, 3))[1];
     218    //Selecting output missing variables with prior information
     219    Matrix _.I = SetCol(_.t);
     220   
    208221    Set _.missingInfo = outputMissingBlock;
     222
     223    //Given initial values of block variables
     224    VMatrix _.store = Mat2VMat(SetMat(Extract(_.missingInfo, 6)));
    209225    Code _.buildWorkSpace = _buildWorkSpace;
    210     VMatrix _.store = Mat2VMat(SetMat(Extract(outputMissingBlock, 6)));
    211226    //Defining the block in the Markov Chain
    212227    Real _defined = define(
    213228      name,
    214       EvalSet(outputMissingBlock, Text(@BSR.MissingBlock b){ b->Name }),
     229      EvalSet(_.missingInfo, Text(@BSR.MissingBlock b){ b->Name }),
    215230      mcmcPartialHandler,
    216231      numBlock,
    217232      firstCol,
    218233      config);
    219     //Output missing pairs i,k
    220     Set _.t = Traspose(Extract(outputMissingBlock, 3))[1];
    221     //Selecting output missing variables with prior information
    222     Matrix _.I = SetCol(_.t);
    223     Set _.prior.normal = Select(outputMissingBlock, Real(@BSR.MissingBlock b)
     234    //Selecting input missing variables with prior information
     235    Set _.prior.normal = Select(_.missingInfo, Real(@BSR.MissingBlock b)
    224236    {
    225237      Or(b->Prior == "Normal",b->Prior == "TruncatedNormal")
    226238    });
    227     //Number of output missing prior equations
     239    //Number of missing prior equations
    228240    Real _.prior.M = Card(_.prior.normal);
    229     //Average vector of output missing priors
    230     Matrix _.prior.nu = If(!_.prior.M, Rand(0,0,0,0),
    231       SetMat(Extract(_.prior.normal, 6)));
    232     //Inverse of standard deviation vector of output missing priors
    233     Matrix _.prior.si = If(!_.prior.M, Rand(0,0,0,0),
    234     {
     241    //Average vector of missing priors
     242    Matrix _.prior.nu = If(!_.prior.M,Rand(0,1,0,0),
     243    {
     244    //WriteLn("TRACE "+_MID+" 3");
     245      SetMat(Extract(_.prior.normal, 6))
     246    });
     247    //Inverse of standard deviation vector of missing priors
     248    Matrix _.prior.si = If(!_.prior.M,Rand(0,1,0,0),
     249    {
     250    //WriteLn("TRACE "+_MID+" 4");
    235251      Matrix s2 = SetMat(Extract(_.prior.normal, 7));
    236252      RPow(s2,-1/2)
    237253    });
    238     //Standarized output of output missing priors
    239     VMatrix _.prior.Y = If(!_.prior.M, Rand(0,1,0,0),
     254    //Standarized output of missing priors
     255    VMatrix _.prior.Y = If(!_.prior.M,Rand(0,1,0,0),
    240256      Mat2VMat(_.prior.nu $* _.prior.si));
    241     //Standarized output of output missing priors
    242     VMatrix _.prior.X = If(!_.prior.M, Rand(0,_.N,0,0),
     257    //Standarized input of missing priors
     258    VMatrix _.prior.X = If(!_.prior.M,Rand(0,_.N,0,0),
    243259    {
    244260    //WriteLn("TRACE "+_MID+" _.prior.X 1");
    … …  
    250266    //WriteLn("TRACE "+_MID+" _.prior.X 3");
    251267      VMatrix pX = Convert(tr,"Cholmod.R.Sparse");
    252     //WriteLn("TRACE "+_MID+" _.prior.X 4 prior.X="<<pX);
     268    //WriteLn("TRACE "+_MID+" _.prior.X 4");
    253269      pX
    254270    })
    255 
    256271  ]])}
    257272
Note: See TracChangeset for help on using the changeset viewer.