Changeset c07b2db in tol


Ignore:
Timestamp:
Apr 1, 2009, 3:52:49 PM (18 years ago)
Author:
Víctor de Buen <vdebuen@…>
Branches:
1941-docker-image-to-test, devel, master
Children:
2dfbf68
Parents:
dd921c2
Message:

Removing max-likelihood step
First approximation to Metropolis-Hastings method explained in #664

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

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

Legend:

Unmodified
Added
Removed

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

    rdd921c2 rc07b2db  
    2828  VMatrix filter(VMatrix X) { X };
    2929  VMatrix draw(Real numSim,NameBlock arima,VMatrix si,VMatrix beta,VMatrix Z)
    30   { _.blockDef::_.store };
    31   Code maxLogLikelihood = draw
     30  { _.blockDef::_.store }
    3231  ]], [[
    3332  //Defining the block in the Markov Chain
    … …  
    4342  //Number of input missing variables (dimension of input missing block)
    4443  Real _.N = Card(_.modelDef->InputMissingBlock);
     44  //Logarithm of density of last generated vector
     45  Real _.betaLogDens = -1/0;
     46
     47  Real _.MH.useAfterIter = Max(1,config::bsr.linBlk.MH.useAfterIter);
     48  Real _.MH.maxRejected  = config::bsr.linBlk.MH.maxRejected;
     49  Real _.MH.maxConsecutive  = config::bsr.linBlk.MH.maxConsecutive;
     50  Real _.MH.tried = 0;
     51  Real _.MH.rejected = 0;
     52  Real _.MH.success = 0;
     53  Real _.MH.consecutive = 0;
     54  Real _.calculatedSimulation = 0;
     55
    4556  //Input missing pairs (equation,variable)
    4657  Matrix _.tj = {
    … …  
    95106    pX
    96107  });
     108
     109  ////////////////////////////////////////////////////////////////////////////
     110  VMatrix getStore(Real unused)
     111  ////////////////////////////////////////////////////////////////////////////
     112  {
     113    _.lin.blk::_.blockDef::_.store
     114  };
     115
    97116  ///////////////////////////////////////////////////////////////////////////
    98117  Real initialize(Real unused)
    … …  
    197216
    198217  ////////////////////////////////////////////////////////////////////////////
     218  Text _.autodoc.member.draw.gibbs =
     219  "Method of Gibss simulation of input missing block";
     220  VMatrix draw.gibbs(Real numSim, NameBlock ls)
     221  ////////////////////////////////////////////////////////////////////////////
     222  {
     223    Real _.MH.consecutive := 0;
     224    Real _.lin.blk::solve(numSim, ls::Y01, ls::X01);
     225    VMatrix _.blockDef::_.store := _.lin.blk::draw(numSim);
     226    VMatrix e = ls::Y01-ls::X01*_.blockDef::_.store;
     227    Real _.betaLogDens := MatSum(Log(f01(VMat2Mat(e))));
     228    _.blockDef::_.store
     229  };
     230
     231  ////////////////////////////////////////////////////////////////////////////
     232  Text _.autodoc.member.draw.metropolis_hastings =
     233  "Method of Metroplois-Hastings simulation of input missing block";
     234  VMatrix draw.metropolis_hastings(Real numSim, NameBlock ls)
     235  ////////////////////////////////////////////////////////////////////////////
     236  {
     237  //WriteLn("TRACE "+_MID+" draw.metropolis_hastings 1");
     238/* * /
     239    WriteLn(" numSim="<<numSim+
     240            " _.MH.success="<<_.MH.success+
     241            " _.MH.tried="<<_.MH.tried+
     242            " _.MH.rejected="<<_.MH.rejected+
     243            " _.MH.consecutive="<<_.MH.consecutive);
     244/* */
     245    Real _.MH.consecutive := _.MH.consecutive + 1;
     246    Real m = VRows(Y);
     247    VMatrix b = getStore(0);
     248    VMatrix ep = ls::Y01-ls::X01*b;
     249    Real p0 = -0.5*VMatDat(MtMSqr(ep),1,1);
     250    Real q0 = _.lin.blk::_.betaLogDens;
     251    Real s0 = Sqrt(-2*p0/m);
     252  //WriteLn(" s0="<<s0+" p0="<<p0+" q0="<<q0);
     253    Real p1 = ?;
     254    Real q1 = ?;
     255    Real accept = False;
     256    Real i = 1;
     257    While(Not(accept) & (i<=_.MH.maxRejected),
     258    {
     259      Real t0 = Copy(Time);
     260      Real _.lin.blk::changed.X := False;
     261      VMatrix b := _.lin.blk::draw(numSim);
     262      VMatrix ep := ls::Y01-ls::X01*b;
     263      Real p1 := -0.5*VMatDat(MtMSqr(ep),1,1);
     264      Real s1 = Sqrt(-2*p1/m);
     265      Real q1 := _.lin.blk::_.betaLogDens;
     266      Real a = Rand(0,1);
     267      Real c = Exp((p1+q0)-(p0+q1));
     268      Real s = Sqrt(gsl_ran_chisq(m)/m);
     269      Real accept := Case(
     270        s1 <= 1, True,
     271        s1 <= s, True,
     272        1  == 1, a<c);
     273      Real tm = Copy(Time) - t0;
     274/* * /
     275      WriteLn(" i="<<i+
     276              " s1="<<s1+
     277              " p1="<<p1+
     278              " q1="<<q1+
     279              " c="<<c+
     280              " a="<<a+
     281              " s="<<s+
     282              " accept="<<accept+
     283              " time="<<tm);
     284/* */
     285      Real _.MH.tried := _.MH.tried +1 ;
     286      Real _.MH.rejected = _.MH.rejected + Not(accept);
     287      Real i := i+1
     288    });
     289    If(accept,
     290    {
     291      Real _.MH.success := _.MH.success + 1;
     292      Real _.betaLogDens := p1;
     293      b
     294    },
     295    {
     296      draw.gibbs(numSim, ls)
     297    }) 
     298  };
     299
     300  ////////////////////////////////////////////////////////////////////////////
    199301  Text _.autodoc.member.draw =
    200302  "Method of generation of input missing block";
    … …  
    202304               NameBlock arima, //ArimaBlock filter
    203305               VMatrix   si,    //SigmaBlock filter
    204                VMatrix   beta,  //Non filtered LinearBlock parameters
    205                VMatrix   Z      //
    206                )
     306               VMatrix   beta,  //Non filtered Main Linear Block parameters
     307               VMatrix   Z)     //Non filtered noise
    207308  ////////////////////////////////////////////////////////////////////////////
    208309  {
    … …  
    210311    {
    211312      NameBlock ls = buildLinearSystem(arima,si,beta,Z);
    212       Real _.lin.blk::solve(numSim, ls::Y01, ls::X01);
    213       VMatrix _.blockDef::_.store := _.lin.blk::draw(numSim);
    214       _.blockDef::_.store
    215     })
    216   };
    217 
    218   ////////////////////////////////////////////////////////////////////////////
    219   Text _.autodoc.member.maxLikelihood =
    220   "Method of contitioned maximum likelihood";
    221   VMatrix maxLikelihood(Real      numSim,
    222                         NameBlock arima, //ArimaBlock filter
    223                         VMatrix   si,    //SigmaBlock filter
    224                         VMatrix   beta,  //Non filtered LinearBlock parameters
    225                         VMatrix   Z      //
    226                )
    227   ////////////////////////////////////////////////////////////////////////////
    228   {
    229     If(Not(_.blockDef::enabled), _.blockDef::_.store,
    230     {
    231       NameBlock ls = buildLinearSystem(arima,si,beta,Z);
    232       VMatrix _.blockDef::_.store:=
    233         _.lin.blk::maxLikelihood(numSim, ls::Y01, ls::X01);
    234       _.blockDef::_.store
     313      Real _.calculatedSimulation := _.calculatedSimulation+1;
     314      If(And(_.calculatedSimulation >_.MH.useAfterIter,
     315             _.MH.consecutive <= _.MH.maxConsecutive),
     316        draw.metropolis_hastings(numSim, ls),
     317        draw.gibbs              (numSim, ls))
    235318    })
    236319  };

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

    rdd921c2 rc07b2db  
    2727  Real initialize(Real unused){ 0 };
    2828  VMatrix filter(VMatrix X) { X };
    29   VMatrix draw(Real numSim,NameBlock arima,VMatrix si,VMatrix beta,VMatrix Z)
    30   { _.blockDef::_.store };
    31   Code maxLogLikelihood = draw
     29  VMatrix draw(Real numSim, VMatrix Y, VMatrix X, NameBlock arm.blk)
     30  { _.blockDef::_.store }
    3231  ]], [[
    3332  //Defining the block in the Markov Chain
    … …  
    4342  //Number of input missing variables (dimension of input missing block)
    4443  Real _.N = Card(varNames);
     44  //Logarithm of density of last generated vector
     45  Real _.betaLogDens = -1/0;
    4546
    46   ///////////////////////////////////////////////////////////////////////////
     47  Real _.MH.useAfterIter = Max(1,config::bsr.linBlk.MH.useAfterIter);
     48  Real _.MH.maxRejected  = config::bsr.linBlk.MH.maxRejected;
     49  Real _.MH.maxConsecutive  = config::bsr.linBlk.MH.maxConsecutive;
     50  Real _.MH.tried = 0;
     51  Real _.MH.rejected = 0;
     52  Real _.MH.success = 0;
     53  Real _.MH.consecutive = 0;
     54  Real _.calculatedSimulation = 0;
     55
     56  ////////////////////////////////////////////////////////////////////////////
    4757  Real initialize(VMatrix beta0, VMatrix A, VMatrix a)
    48   ///////////////////////////////////////////////////////////////////////////
     58  ////////////////////////////////////////////////////////////////////////////
    4959  {
     60   
    5061    _.lin.blk::initialize(beta0, A, a)
    5162  };
    … …  
    6071
    6172  ////////////////////////////////////////////////////////////////////////////
    62   Text _.autodoc.member.draw =
    63   "Method of generation of mian linear block";
    64   VMatrix draw(Real numSim, VMatrix Y, VMatrix X)
     73  Text _.autodoc.member.draw.gibbs =
     74  "Method of Gibbs simulation of main linear block";
     75  VMatrix draw.gibbs(Real numSim, VMatrix Y, VMatrix X, NameBlock arm.blk)
    6576  ////////////////////////////////////////////////////////////////////////////
    6677  {
    67     Real _.lin.blk::solve(numSim, Y, X);
    68     _.lin.blk::draw(numSim)
     78  //WriteLn("TRACE "+_MID+" draw.gibbs 1");
     79    Real _.MH.consecutive := 0;
     80    VMatrix Y01 = arm.blk::filter(Y);
     81    VMatrix X01 = arm.blk::filter(X);
     82    Real _.lin.blk::solve(numSim, Y01, X01);
     83    VMatrix b =_.lin.blk::draw(numSim);
     84    VMatrix ep = Y01-X01*b;
     85    Real _.betaLogDens := -0.5*VMatDat(MtMSqr(ep),1,1);
     86    b
    6987  };
    7088
    7189  ////////////////////////////////////////////////////////////////////////////
    72   Text _.autodoc.member.maxLikelihood =
    73   "Method of contitioned maximum likelihood";
    74   VMatrix maxLikelihood(Real      numSim,
    75                         NameBlock arima, //ArimaBlock filter
    76                         VMatrix   si,    //SigmaBlock filter
    77                         VMatrix   beta,  //Non filtered LinearBlock parameters
    78                         VMatrix   Z      //
    79                )
     90  Text _.autodoc.member.draw.metropolis_hastings =
     91  "Method of Metroplois-Hastings simulation of main linear block";
     92  VMatrix draw.metropolis_hastings(Real numSim, VMatrix Y, VMatrix X, NameBlock arm.blk)
     93  ////////////////////////////////////////////////////////////////////////////
     94  {
     95  //WriteLn("TRACE "+_MID+" draw.metropolis_hastings 1");
     96/* * /
     97    WriteLn(" numSim="<<numSim+
     98            " _.MH.success="<<_.MH.success+
     99            " _.MH.tried="<<_.MH.tried+
     100            " _.MH.rejected="<<_.MH.rejected+
     101            " _.MH.consecutive="<<_.MH.consecutive);
     102/* */
     103    Real _.MH.consecutive := _.MH.consecutive + 1;
     104    Real m = VRows(Y);
     105    VMatrix b = getStore(0);
     106    VMatrix ep = arm.blk::filter(Y-X*b);
     107    Real p0 = -0.5*VMatDat(MtMSqr(ep),1,1);
     108    Real q0 = _.lin.blk::_.betaLogDens;
     109    Real s0 = Sqrt(-2*p0/m);
     110  //WriteLn(" s0="<<s0+" p0="<<p0+" q0="<<q0);
     111    Real p1 = ?;
     112    Real q1 = ?;
     113    Real accept = False;
     114    Real i = 1;
     115    While(Not(accept) & (i<=_.MH.maxRejected),
     116    {
     117      Real t0 = Copy(Time);
     118      Real _.lin.blk::changed.X := False;
     119      VMatrix b := _.lin.blk::draw(numSim);
     120      VMatrix ep := arm.blk::filter(Y-X*b);
     121      Real p1 := -0.5*VMatDat(MtMSqr(ep),1,1);
     122      Real s1 = Sqrt(-2*p1/m);
     123      Real q1 := _.lin.blk::_.betaLogDens;
     124      Real a = Rand(0,1);
     125      Real c = Exp((p1+q0)-(p0+q1));
     126      Real s = Sqrt(gsl_ran_chisq(m)/m);
     127      Real accept := Case(
     128        s1 <= 1, True,
     129        s1 <= s, True,
     130        1  == 1, a<c);
     131      Real tm = Copy(Time) - t0;
     132/* * /
     133      WriteLn(" i="<<i+
     134              " s1="<<s1+
     135              " p1="<<p1+
     136              " q1="<<q1+
     137              " c="<<c+
     138              " a="<<a+
     139              " s="<<s+
     140              " accept="<<accept+
     141              " time="<<tm);
     142/* */
     143      Real _.MH.tried := _.MH.tried +1 ;
     144      Real _.MH.rejected = _.MH.rejected + Not(accept);
     145      Real i := i+1
     146    });
     147    If(accept,
     148    {
     149      Real _.MH.success := _.MH.success + 1;
     150      Real _.betaLogDens := p1;
     151      b
     152    },
     153    {
     154      draw.gibbs(numSim, Y, X, arm.blk)
     155    }) 
     156  };
     157
     158  ////////////////////////////////////////////////////////////////////////////
     159  Text _.autodoc.member.draw =
     160  "Method of generation of mian linear block";
     161  VMatrix draw(Real numSim, VMatrix Y, VMatrix X, NameBlock arm.blk)
    80162  ////////////////////////////////////////////////////////////////////////////
    81163  {
    82164    If(Not(_.blockDef::enabled), _.blockDef::_.store,
    83165    {
    84       NameBlock ls = buildLinearSystem(arima,si,beta,Z);
    85       VMatrix _.blockDef::_.store:=
    86         _.lin.blk::maxLikelihood(numSim, ls::Y01, ls::X01);
    87       _.blockDef::_.store
    88     })
     166      Real _.calculatedSimulation := _.calculatedSimulation+1;
     167      If(And(_.calculatedSimulation >_.MH.useAfterIter,
     168             _.MH.consecutive <= _.MH.maxConsecutive),
     169        draw.metropolis_hastings(numSim, Y, X, arm.blk),
     170        draw.gibbs              (numSim, Y, X, arm.blk))
     171   })
    89172  };
    90173

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

    rdd921c2 rc07b2db  
    2727  VMatrix filter(VMatrix Y) { Y };
    2828  VMatrix draw(Real numSim, NameBlock arima, VMatrix si, VMatrix Z)
    29   { _.blockDef::_.store };
    30   Code maxLogLikelihood = draw
     29  { _.blockDef::_.store }
    3130  ]],[[
    3231  //Defining the block in the Markov Chain
    … …  
    4241  //Number of output missing variables (dimension of output missing block)
    4342  Real _.N = Card(_.modelDef->OutputMissingBlock);
     43  //Logarithm of density of last generated vector
     44  Real _.betaLogDens = -1/0;
     45
     46  Real _.MH.useAfterIter = Max(1,config::bsr.linBlk.MH.useAfterIter);
     47  Real _.MH.maxRejected  = config::bsr.linBlk.MH.maxRejected;
     48  Real _.MH.maxConsecutive  = config::bsr.linBlk.MH.maxConsecutive;
     49  Real _.MH.tried = 0;
     50  Real _.MH.rejected = 0;
     51  Real _.MH.success = 0;
     52  Real _.MH.consecutive = 0;
     53  Real _.calculatedSimulation = 0;
     54
    4455  //Output missing pairs i,k
    4556  Set _.t = Traspose(Extract(outputMissingBlock, 3))[1];
    … …  
    7889    pX
    7990  });
     91
     92  ////////////////////////////////////////////////////////////////////////////
     93  VMatrix getStore(Real unused)
     94  ////////////////////////////////////////////////////////////////////////////
     95  {
     96    _.lin.blk::_.blockDef::_.store
     97  };
     98
    8099  ///////////////////////////////////////////////////////////////////////////
    81100  Real initialize(Real unused)
    … …  
    179198
    180199  ////////////////////////////////////////////////////////////////////////////
     200  Text _.autodoc.member.draw.gibbs =
     201  "Method of Gibss simulation of output missing block";
     202  VMatrix draw.gibbs(Real numSim, NameBlock ls)
     203  ////////////////////////////////////////////////////////////////////////////
     204  {
     205    Real _.MH.consecutive := 0;
     206    Real _.lin.blk::solve(numSim, ls::Y01, ls::X01);
     207    VMatrix _.blockDef::_.store := _.lin.blk::draw(numSim);
     208    VMatrix e = ls::Y01-ls::X01*_.blockDef::_.store;
     209    Real _.betaLogDens := MatSum(Log(f01(VMat2Mat(e))));
     210    _.blockDef::_.store
     211  };
     212
     213  ////////////////////////////////////////////////////////////////////////////
     214  Text _.autodoc.member.draw.metropolis_hastings =
     215  "Method of Metroplois-Hastings simulation of output missing block";
     216  VMatrix draw.metropolis_hastings(Real numSim, NameBlock ls)
     217  ////////////////////////////////////////////////////////////////////////////
     218  {
     219  //WriteLn("TRACE "+_MID+" draw.metropolis_hastings 1");
     220/* * /
     221    WriteLn(" numSim="<<numSim+
     222            " _.MH.success="<<_.MH.success+
     223            " _.MH.tried="<<_.MH.tried+
     224            " _.MH.rejected="<<_.MH.rejected+
     225            " _.MH.consecutive="<<_.MH.consecutive);
     226/* */
     227    Real _.MH.consecutive := _.MH.consecutive + 1;
     228    Real m = VRows(Y);
     229    VMatrix b = getStore(0);
     230    VMatrix ep = ls::Y01-ls::X01*b;
     231    Real p0 = -0.5*VMatDat(MtMSqr(ep),1,1);
     232    Real q0 = _.lin.blk::_.betaLogDens;
     233    Real s0 = Sqrt(-2*p0/m);
     234  //WriteLn(" s0="<<s0+" p0="<<p0+" q0="<<q0);
     235    Real p1 = ?;
     236    Real q1 = ?;
     237    Real accept = False;
     238    Real i = 1;
     239    While(Not(accept) & (i<=_.MH.maxRejected),
     240    {
     241      Real t0 = Copy(Time);
     242      Real _.lin.blk::changed.X := False;
     243      VMatrix b := _.lin.blk::draw(numSim);
     244      VMatrix ep := ls::Y01-ls::X01*b;
     245      Real p1 := -0.5*VMatDat(MtMSqr(ep),1,1);
     246      Real s1 = Sqrt(-2*p1/m);
     247      Real q1 := _.lin.blk::_.betaLogDens;
     248      Real a = Rand(0,1);
     249      Real c = Exp((p1+q0)-(p0+q1));
     250      Real s = Sqrt(gsl_ran_chisq(m)/m);
     251      Real accept := Case(
     252        s1 <= 1, True,
     253        s1 <= s, True,
     254        1  == 1, a<c);
     255      Real tm = Copy(Time) - t0;
     256/* * /
     257      WriteLn(" i="<<i+
     258              " s1="<<s1+
     259              " p1="<<p1+
     260              " q1="<<q1+
     261              " c="<<c+
     262              " a="<<a+
     263              " s="<<s+
     264              " accept="<<accept+
     265              " time="<<tm);
     266/* */
     267      Real _.MH.tried := _.MH.tried +1 ;
     268      Real _.MH.rejected = _.MH.rejected + Not(accept);
     269      Real i := i+1
     270    });
     271    If(accept,
     272    {
     273      Real _.MH.success := _.MH.success + 1;
     274      Real _.betaLogDens := p1;
     275      b
     276    },
     277    {
     278      draw.gibbs(numSim, ls)
     279    }) 
     280  };
     281
     282  ////////////////////////////////////////////////////////////////////////////
    181283  Text _.autodoc.member.draw =
    182284  "Method of generation of output missing block";
    … …  
    191293    {
    192294      NameBlock ls = buildLinearSystem(arima,si,Z);
    193       Real _.lin.blk::solve(numSim, ls::Y01, ls::X01);
    194       VMatrix _.blockDef::_.store := _.lin.blk::draw(numSim);
    195       _.blockDef::_.store
     295      Real _.calculatedSimulation := _.calculatedSimulation+1;
     296      If(And(_.calculatedSimulation >_.MH.useAfterIter,
     297             _.MH.consecutive <= _.MH.maxConsecutive),
     298        draw.metropolis_hastings(numSim, ls),
     299        draw.gibbs              (numSim, ls))
    196300    })
    197301  };
    198 
    199   ////////////////////////////////////////////////////////////////////////////
    200   Text _.autodoc.member.maxLikelihood =
    201   "Method of contitioned maximum likelihood";
    202   VMatrix maxLikelihood(Real      numSim,
    203                         NameBlock arima, //ArimaBlock filter
    204                         VMatrix   si,    //SigmaBlock filter
    205                         VMatrix   Z      //
    206                )
    207   ////////////////////////////////////////////////////////////////////////////
    208   {
    209     If(Not(_.blockDef::enabled), _.blockDef::_.store,
    210     {
    211       NameBlock ls = buildLinearSystem(arima,si,Z);
    212       VMatrix _.blockDef::_.store :=
    213         _.lin.blk::maxLikelihood(numSim, ls::Y01, ls::X01);
    214       _.blockDef::_.store
    215     })
    216   };
    217 
    218302
    219303  ////////////////////////////////////////////////////////////////////////////
Note: See TracChangeset for help on using the changeset viewer.