Changeset 74279bb3 in tol


Ignore:
Timestamp:
Sep 4, 2009, 5:35:09 PM (17 years ago)
Author:
Víctor de Buen <vdebuen@…>
Branches:
1941-docker-image-to-test, devel, master
Children:
b8a0836
Parents:
9efcb47
Message:

Implementing partial simulation for sigma block

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

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r9efcb47 r74279bb3  
    33// MODULE: Stat (Statistics)
    44// FILE: _sigmaBlock.tol
    5 // PURPOSE: Declares method BysMcmc::Bsr::Gibbs::SigmaBlock
     5// PURPOSE: Declares Class BysMcmc::Bsr::Gibbs::BlockSigma
    66/////////////////////////////////////////////////////////////////////////////
     7
     8//--------------------------------------------------------------------------
     9Struct BlockSigma.PartialIndex
     10//--------------------------------------------------------------------------
     11{
     12  Real idxFixed;
     13  Real idxReload
     14};
    715
    816/////////////////////////////////////////////////////////////////////////////
    … …  
    3139  //--------------------------------------------------------------------------
    3240  NameBlock _.ws = [[ Real unused_=? ]];
     41  Set _.segment.partial.index = Copy(Empty);
    3342
    3443  //--------------------------------------------------------------------------
    3544  //Methods
    3645  //--------------------------------------------------------------------------
     46
     47  ////////////////////////////////////////////////////////////////////////////
     48  Real initialize(Real unused)
     49  ////////////////////////////////////////////////////////////////////////////
     50  {
     51  //WriteLn("TRACE "+_MID+" initialize 1 ");
     52    Real hasFixed   = _.partialSampler::_.fixedParamNum  > 0;
     53  //WriteLn("TRACE "+_MID+" initialize hasFixed="<<hasFixed);
     54    Real hasReload  = _.partialSampler::_.reloadParamNum > 0;
     55  //WriteLn("TRACE "+_MID+" initialize hasReload="<<hasReload);
     56    Real hasPartial = hasFixed | hasReload;
     57  //WriteLn("TRACE "+_MID+" initialize hasPartial="<<hasPartial);
     58    Set If(hasPartial,
     59    {
     60      Set _.segment.partial.index :=
     61      EvalSet(_.sigmaBlockVar,Set(BSR.NoiseDistrib segment)
     62      {
     63        Real idxFixed  = 0;
     64        Real idxReload = 0;
     65        Text name = ToName(segment->SigmaName);
     66        Real If(hasFixed,
     67        {
     68          idxFixed:=FindIndexByName(_.partialSampler::_.fixedParam, name)
     69        });
     70        Real If(hasReload & !idxFixed,
     71        {
     72          idxReload:=FindIndexByName(_.partialSampler::_.reloadParam, name)
     73        });
     74        BysMcmc::Bsr::Gibbs::BlockSigma.PartialIndex(idxFixed, idxReload)
     75      })
     76    });
     77    1
     78  };
     79
    3780  ////////////////////////////////////////////////////////////////////////////
    3881  Real _buildWorkSpace(VMatrix residuals_)
    … …  
    78121
    79122  ////////////////////////////////////////////////////////////////////////////
     123  Real draw.invRandChisq.segment(
     124    Real numSim,
     125    VMatrix residuals,
     126    BSR.NoiseDistrib segment)
     127  ////////////////////////////////////////////////////////////////////////////
     128  {
     129  //WriteLn("TRACE "+_MID+" draw.invRandChisq 2.1 "+segment->Name);
     130    NameBlock sumSqr = sumOfSquares(numSim, residuals, segment);
     131  //WriteLn("TRACE "+_MID+" draw.invRandChisq 2.2 ");
     132    Case(
     133    StructName(segment->SigmaPrior)=="Bsr.Sigma.Prior.InverseGamma",
     134    {
     135    //WriteLn("TRACE "+_MID+" draw.invRandChisq 2.2.1 ");
     136      Real s2 = segment->SigmaPrior->Sigma2;
     137    //WriteLn("TRACE "+_MID+" draw.invRandChisq 2.2.2 ");
     138      Real w  = segment->SigmaPrior->Weight;
     139    //WriteLn("TRACE ["+_MID+"] draw.invRandChisq("<<s2+","<<w+")");
     140      Case(
     141      w==1, 
     142      {
     143      //WriteLn("TRACE "+_MID+" draw.invRandChisq 2.2.3.1 ");
     144        s2
     145      },
     146      w==0,
     147      {
     148      //WriteLn("TRACE "+_MID+" draw.invRandChisq 2.2.3.2 ");
     149        sumSqr::rss/RandChisq(sumSqr::m)
     150      },
     151      LT(0,w,1),
     152      {
     153      //WriteLn("TRACE "+_MID+" draw.invRandChisq 2.2.3.3 ");
     154        //w = m0/(m0+m); m0*(w-1)=-m*w; m0 = m*w/(1-w)
     155        Real m0 = Max(1,Round(sumSqr::m*w/(1-w)));
     156      //WriteLn("TRACE "+_MID+" draw.invRandChisq 2.2.3.4 ");
     157        Real rss = sumSqr::rss + m0 * s2;
     158      //WriteLn("TRACE "+_MID+" draw.invRandChisq 2.2.3.5 ");
     159        Real m = sumSqr::m + m0;
     160      //WriteLn("TRACE "+_MID+" draw.invRandChisq 2.2.3.6 ");
     161        Real S2 = rss/RandChisq(m);
     162      //WriteLn("TRACE "+_MID+" draw.invRandChisq 2.2.3.7 ");
     163        /*
     164      //WriteLn("TRACE "+_MID+" draw.invRandChisq ("+segment->Name+") "+
     165          " s2="<<s2+
     166          " w="<<w+
     167          " m0="<<m0+
     168          " sumSqr::m="<<sumSqr::m+
     169          " sumSqr::rss="<<sumSqr::rss+
     170          " m="<<m+
     171          " S2="<<S2); */
     172        S2         
     173      },
     174      1==1,
     175      {
     176        WriteLn(_MID+" Weigtht of prior sigma is not in [0,1] "
     177        "for segment "+segment->Name,"E");
     178        ?
     179      })
     180
     181    },
     182    1==1,
     183    {
     184    //WriteLn("TRACE "+_MID+" draw.invRandChisq 2.3 ");
     185      sumSqr::rss/RandChisq(sumSqr::m)
     186    })
     187  };
     188
     189  ////////////////////////////////////////////////////////////////////////////
    80190  Text _.autodoc.member.draw.invRandChisq =
    81191  "Method of generation of sigma block containing draws of all variant "
    … …  
    85195  {
    86196  //WriteLn("TRACE "+_MID+" draw.invRandChisq 1 ");
    87     Set v = EvalSet(_.sigmaBlockVar,Real(BSR.NoiseDistrib segment)
    88     {
    89     //WriteLn("TRACE "+_MID+" draw.invRandChisq 2.1 "+segment->Name);
    90       NameBlock sumSqr = sumOfSquares(numSim, residuals, segment);
    91     //WriteLn("TRACE "+_MID+" draw.invRandChisq 2.2 ");
    92       Case(
    93       StructName(segment->SigmaPrior)=="Bsr.Sigma.Prior.InverseGamma",
    94       {
    95       //WriteLn("TRACE "+_MID+" draw.invRandChisq 2.2.1 ");
    96         Real s2 = segment->SigmaPrior->Sigma2;
    97       //WriteLn("TRACE "+_MID+" draw.invRandChisq 2.2.2 ");
    98         Real w  = segment->SigmaPrior->Weight;
    99       //WriteLn("TRACE ["+_MID+"] draw.invRandChisq("<<s2+","<<w+")");
    100         Case(
    101         w==1, 
    102         {
    103         //WriteLn("TRACE "+_MID+" draw.invRandChisq 2.2.3.1 ");
    104           s2
     197   
     198    Set v = If(!Card(_.segment.partial.index),
     199    {
     200      EvalSet(_.sigmaBlockVar,Real(BSR.NoiseDistrib segment)
     201      {
     202        draw.invRandChisq.segment(numSim, residuals, segment)
     203      })
     204    },
     205    {
     206      For(1,Card(_.sigmaBlockVar),Real(Real k)
     207      {
     208        Set partial = _.segment.partial.index[k];
     209        Real Case(partial->idxFixed,
     210        {
     211          VMatDat(_.partialSampler::_.fixedValues,partial->idxFixed,1)
    105212        },
    106         w==0,
    107         {
    108         //WriteLn("TRACE "+_MID+" draw.invRandChisq 2.2.3.2 ");
    109           sumSqr::rss/RandChisq(sumSqr::m)
     213        {
     214          VMatDat(_.partialSampler::_.reloadValues,partial->idxreload,1)
    110215        },
    111         LT(0,w,1),
    112         {
    113         //WriteLn("TRACE "+_MID+" draw.invRandChisq 2.2.3.3 ");
    114           //w = m0/(m0+m); m0*(w-1)=-m*w; m0 = m*w/(1-w)
    115           Real m0 = Max(1,Round(sumSqr::m*w/(1-w)));
    116         //WriteLn("TRACE "+_MID+" draw.invRandChisq 2.2.3.4 ");
    117           Real rss = sumSqr::rss + m0 * s2;
    118         //WriteLn("TRACE "+_MID+" draw.invRandChisq 2.2.3.5 ");
    119           Real m = sumSqr::m + m0;
    120         //WriteLn("TRACE "+_MID+" draw.invRandChisq 2.2.3.6 ");
    121           Real S2 = rss/RandChisq(m);
    122         //WriteLn("TRACE "+_MID+" draw.invRandChisq 2.2.3.7 ");
    123           /*
    124         //WriteLn("TRACE "+_MID+" draw.invRandChisq ("+segment->Name+") "+
    125             " s2="<<s2+
    126             " w="<<w+
    127             " m0="<<m0+
    128             " sumSqr::m="<<sumSqr::m+
    129             " sumSqr::rss="<<sumSqr::rss+
    130             " m="<<m+
    131             " S2="<<S2); */
    132           S2         
    133         },
    134         1==1,
    135         {
    136           WriteLn(_MID+" Weigtht of prior sigma is not in [0,1] "
    137           "for segment "+segment->Name,"E");
    138           ?
     216        {
     217          draw.invRandChisq.segment(numSim, residuals, _.sigmaBlockVar[k])
    139218        })
    140 
    141       },
    142       1==1,
    143       {
    144       //WriteLn("TRACE "+_MID+" draw.invRandChisq 2.3 ");
    145         sumSqr::rss/RandChisq(sumSqr::m)
    146219      })
    147220    });
Note: See TracChangeset for help on using the changeset viewer.