Changeset f1fd2aa in tol


Ignore:
Timestamp:
May 13, 2010, 7:17:07 AM (16 years ago)
Author:
Víctor de Buen <vdebuen@…>
Branches:
1941-docker-image-to-test, devel, master
Children:
97a3b1b
Parents:
6d9639b
Message:

El método accept-reject se activa cuando hay restricciones de desigualdad arbitrarias y está controlado por el miembro maxAcceptRejectTry

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

File:
1 edited

Legend:

Unmodified
Added
Removed

  • tol_pkg/BysVecLinReg/Multivariate.tol

    r6d9639b rf1fd2aa  
    6262  //Stacked columns of A (D*Mx1)
    6363  VMatrix _.vec.A = Rand(0,0,0,0);
     64  VMatrix _.vec.A.prior.nu = Rand(0,0,0,0);
     65  VMatrix _.vec.A.prior.S = Rand(0,0,0,0);
    6466  //Average of posterior distribution of A for alpha=0: Y*X'*(X*Xt)^-1 (dxm)
    6567  VMatrix _.nu = Rand(0,0,0,0);
    … …  
    8789  VMatrix _.VxIm.Le = Rand(0,0,0,0);
    8890
     91  Real maxAcceptRejectTry = 100;
     92  Real _.totalTry = 0;
     93
    8994  ////////////////////////////////////////////////////////////////////////////
    9095  // Optional linear constraining equalities
    … …  
    98103  //Right vector of linear constraining equalities
    99104  VMatrix _.c = Rand(0,0,0,0);
     105  //Constraining equalities in standariced metric used internally
     106  VMatrix _.C01 = Rand(0,0,0,0);
     107  VMatrix _.c01 = Rand(0,0,0,0);
    100108  ////////////////////////////////////////////////////////////////////////////
    101109
    … …  
    143151      VMatrix _.C := C;
    144152      VMatrix _.c := c;
    145       True
     153      recalc.C01c01(0)
    146154    })
    147155  };
    … …  
    151159  ////////////////////////////////////////////////////////////////////////////
    152160  {
     161    Real _.totalTry := 0;
    153162    Real _.d := VRows(Y);
    154163    Real _.m := VRows(X);
    … …  
    172181    1==1,
    173182    {
    174       change.X(X);
    175       change.Y(Y);
    176       change.alpha(0);
     183      Real change.X(X);
     184      Real change.Y(Y);
     185      Real change.alpha(0);
    177186      VMatrix _.vec.A := (_.vec.nu)/(1+_.alpha);
    178       VMatrix _.A := vec2vmatrix(_.vec.A,Constant(_.d,_.m,1));
    179       change.V(Eye(_.d));
     187      VMatrix _.A := BysVecLinReg::vec2vmatrix(_.vec.A,Constant(_.d,_.m,1));
     188      Real change.V(Eye(_.d));
     189      VMatrix _.vec.A.prior.nu := _.vec.nu*0;
     190      VMatrix _.vec.A.prior.S := _.V*0;
    180191      True
    181192    })
    … …  
    241252      VMatrix _.nu:=Tra(CholeskiSolve(_.X.Lf,XYt,"PtLLtP"));
    242253      VMatrix _.nuXYt := _.nu*XYt;
    243       VMatrix _.vec.nu:=vmatrix2vec(_.nu);
     254      VMatrix _.vec.nu:=BysVecLinReg::vmatrix2vec(_.nu);
    244255      True
    245256    })
    … …  
    263274
    264275  ////////////////////////////////////////////////////////////////////////////
     276  Real recalc.C01c01(Real void)
     277  ////////////////////////////////////////////////////////////////////////////
     278  {
     279    VMatrix _.C01 := Tra(_.IdxXLi*Tra(_.C))*_.VxIm.Le;
     280    VMatrix _.c01 := _.c - _.C * _.vec.nu;
     281    True
     282  };
     283
     284  ////////////////////////////////////////////////////////////////////////////
    265285  Real change.V(VMatrix V)
    266286  ////////////////////////////////////////////////////////////////////////////
    … …  
    269289    {
    270290      VMatrix _.V := (V+Tra(V))/2;
     291    //WriteLn("TRACE change.V 1");
    271292      VMatrix _.VxIm := KroneckerProduct(_.V, Eye(_.m));
    272293      VMatrix _.VxIm.Lf := CholeskiFactor(_.VxIm,"X");
    … …  
    276297      VMatrix _.VxIm.Le := CholeskiSolve(_.VxIm.Lf,
    277298        Convert(_.VxIm.Lf, "Cholmod.R.Sparse"),"P");
    278       True
    279     })
    280   };
     299      If(!_.r.C, True, recalc.C01c01(0))
     300    })
     301  };
     302
    281303
    282304  ////////////////////////////////////////////////////////////////////////////
    … …  
    309331  {
    310332    VMatrix gibbs.draw.V(numSim);
    311     vmatrix2vec(_.V)
     333    BysVecLinReg::vmatrix2vec(_.V)
    312334  };
    313335
    … …  
    316338  ////////////////////////////////////////////////////////////////////////////
    317339  {
    318     If(!_.r.C, Gaussian(_.d*_.m,1,0,1),
    319     {
    320       VMatrix C = Tra(_.IdxXLi*Tra(_.C))*_.VxIm.Le;
    321       VMatrix c = _.c - _.C * _.vec.nu;
    322       eqConstrainedStdGaussian(C,c)
    323     })
    324   };
    325 
    326   ////////////////////////////////////////////////////////////////////////////
    327   VMatrix gibbs.draw.vec.A(Real numSim)
    328   ////////////////////////////////////////////////////////////////////////////
    329   {
     340    If(!_.r.C,
     341      Gaussian(_.d*_.m,1,0,1),
     342      BysVecLinReg::eqConstrainedStdGaussian(_.C01,_.c01))
     343  };
     344
     345  ////////////////////////////////////////////////////////////////////////////
     346  VMatrix try.draw.vec.A(Real numSim)
     347  ////////////////////////////////////////////////////////////////////////////
     348  {
     349    Real _.totalTry := _.totalTry + 1;
    330350    VMatrix e = gibbs.draw.vec.A.std(numSim);
    331351    VMatrix e.1 = _.VxIm.Le*e;
    332352    VMatrix e.2 = Tra(Tra(e.1)*_.IdxXLi);
    333     VMatrix _.vec.A := (_.vec.nu + e.2)/(1+_.alpha);
    334     _.vec.A
     353    (_.vec.nu + e.2)/(1+_.alpha)
     354  };
     355
     356  ////////////////////////////////////////////////////////////////////////////
     357  VMatrix gibbs.draw.vec.A(Real numSim)
     358  ////////////////////////////////////////////////////////////////////////////
     359  {
     360  //WriteLn("TRACE gibbs.draw.vec.A 1");
     361    VMatrix vec.A = try.draw.vec.A(numSim);
     362    If(!_.r.G,
     363    {
     364    //WriteLn("TRACE gibbs.draw.vec.A 2");
     365      VMatrix _.vec.A := vec.A;
     366      vec.A
     367    },
     368    {
     369      Real match = False;
     370      Real _numTry = 0;
     371      VMatrix A = BysVecLinReg::vec2vmatrix(vec.A,_.A);
     372      While(And(!match,_numTry<maxAcceptRejectTry),
     373      {
     374      //WriteLn("TRACE gibbs.draw.vec.A 3");
     375        VMatrix G.A = eval.arbitraryConstrainingInequalities(A);
     376        Real match := VMatMin(G.A)>=0;
     377        VMatrix If(!match,
     378        {
     379          VMatrix vec.A := try.draw.vec.A(numSim);
     380          VMatrix A     := BysVecLinReg::vec2vmatrix(vec.A,_.A)
     381        });
     382        Real _numTry := _numTry+1
     383      });
     384      If(match,
     385      {
     386        VMatrix _.vec.A := vec.A;
     387        VMatrix _.A := A;
     388        Copy(_.vec.A)
     389      },
     390      {
     391        WriteLn("[BysVecLinReg::@Multivariate::gibbs.draw.vec.A] Too much "
     392        "rejected tries at simulation "<<numSim, "W");
     393        Copy(_.vec.A)
     394      })
     395    })
    335396  };
    336397
    … …  
    340401  {
    341402    VMatrix gibbs.draw.vec.A(numSim);
    342     VMatrix _.A := vec2vmatrix(_.vec.A,_.A);
     403    If(!_.r.G,
     404    {
     405      VMatrix _.A := BysVecLinReg::vec2vmatrix(_.vec.A,_.A)
     406    });
    343407    Copy(_.A)
    344408  };
    … …  
    348412  ////////////////////////////////////////////////////////////////////////////
    349413  {
    350     vmatrix2vec.getPosition(i,j,_.d,_.m)
     414    BysVecLinReg::vmatrix2vec.getPosition(i,j,_.d,_.m)
    351415  }
    352416
Note: See TracChangeset for help on using the changeset viewer.