Changeset dcc06082 in tol for tol_tests


Ignore:
Timestamp:
Apr 20, 2010, 2:17:43 PM (17 years ago)
Author:
Víctor de Buen <vdebuen@…>
Branches:
1941-docker-image-to-test, devel, master
Children:
81a5a9b
Parents:
53196af
Message:

massive tests

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

Location:
tol_tests/tol/ARIMA
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • tol_tests/tol/ARIMA/Stationarity/test_rand_stationary.tol

    r53196af rdcc06082  
    22Matrix rs2 = Group("ConcatRows",For(1,1000,Matrix(Real n)
    33{
    4   -SubCol(Tra(PolMat(RandStationary(2),3,1)),[[2,3]])
     4  SubCol(Tra(PolMat(RandStationary(2),3,1)),[[2,3]])
    55}));

  • tol_tests/tol/ARIMA/Stationarity/test_uniform_rand_stationary.tol

    r53196af rdcc06082  
    11
    2 #Require MatQuery;
    3 
    4 
    5 //////////////////////////////////////////////////////////////////////////////
    6 Set getFaces(Matrix sample)
    7 //////////////////////////////////////////////////////////////////////////////
     2Real maxIterByTry = 20;
     3Real deg = 3;
     4Real size = 10000;
     5Real total_success = 0;
     6Real total_tries = 0;
     7Real total_iter = 0;
     8Matrix sample = Rand(0,deg,0,0);
     9Set rng_col = Range(2,deg+1,1);
     10While(total_success<size,
    811{
    9   Real c = Columns(sample);
    10   SetConcat(For(1,c-1,Set(Real i)
     12  Polyn a = RandStationary(deg);
     13  Polyn b = RandStationary(deg);
     14//Matrix v = Gaussian(1,deg,0,1);
     15  Matrix v = SubCol(Tra(PolMat(a-b,deg+1,1)),rng_col);
     16  Real norm = Sqrt(MatDat(MtMSqr(v),1,1));
     17  Matrix u = v * (1/norm);
     18  Polyn U = MatPol(u)*B;
     19  Polyn p = Copy(a);
     20  Real iter = 0;
     21  Real ok = False;
     22  Real total_tries := total_tries+1;
     23  While(And(!ok,iter<maxIterByTry),
    1124  {
    12     For(i+1,c,Matrix(Real j)
    13     {
    14       Matrix face = SubCol(sample,[[i,j]]);
    15       Eval("Face_"<<i+"_"<<j+"=face")
    16     })
    17   }))
    18 };
    19 
    20 Real maxIterByTry = 0;
    21 Real n = 3;
    22 Real size = 2000;
    23 Real distance.target = 0.001;
    24 Real neighbour.num = 10;
    25 
    26 Real candidate.factor = 1;
    27 
    28 //////////////////////////////////////////////////////////////////////////////
    29 Matrix draw.arbitrary(Real void)
    30 //////////////////////////////////////////////////////////////////////////////
    31 {
    32   -SubCol(Tra(PolMat(RandStationary(n),n+1,1)),Range(2,n+1,1))
    33 };
    34 
    35 
    36 //////////////////////////////////////////////////////////////////////////////
    37 Real isInDomain(Matrix x)
    38 //////////////////////////////////////////////////////////////////////////////
    39 {
    40   IsStationary(1-MatPol(x)*B)
    41 };
    42 
    43 Matrix candidates =
    44 {
    45   Matrix sample = Group("ConcatRows",For(1,size,draw.arbitrary));
    46   MatQuery::@ANN.KDTree ann = MatQuery::@ANN.KDTree::New(sample);
    47   MatQuery::@KNN.Neighbourhood neighbour = ann::search(sample,2);
    48   Real distance.average = MatAvr(Sqrt(SubCol(neighbour::_.distance2,[[2]])));
    49   Real candidate.factor := Ceil(distance.average/distance.target);
    50   Matrix If(candidate.factor>1,sample,
    51     sample := sample <<
    52     Group("ConcatRows",For(1,size*(candidate.factor-1),draw.arbitrary)));
    53   sample << Group("ConcatRows",For(2,neighbour.num,Matrix(Real k)
    54   {
    55     sample + Rand(size,n,-distance.target,+distance.target)
    56   }))
    57 };
    58 
    59 Real candidate.num = Rows(candidates);
    60 
    61 Matrix candidates.isInDomain = SetCol(For(1,candidate.num,Real(Real i)
    62 {
    63   isInDomain(SubRow(candidates,[[i]]))
    64 }));
    65 
    66 Real candidates.isInDomain.ratio = MatAvr(candidates.isInDomain);
    67 
    68 Real numOutDomain = candidate.num-MatSum(candidates.isInDomain);
    69 
    70 
    71 Matrix candidates.sample = SubRow(candidates,
    72   MatSet(Sub(RandPermutation(1,Rows(candidates)),1,1,1,size))[1]);
    73 
    74 Set candidates.showFaces = getFaces(candidates.sample);
    75 
    76 MatQuery::@ANN.KDTree ann = MatQuery::@ANN.KDTree::New(candidates);
    77 
    78 MatQuery::@KNN.Neighbourhood neighbour = ann::search(candidates,neighbour.num);
    79 
    80 Real distance.average = MatAvr(Sqrt(SubCol(neighbour::_.distance2,[[2]])));
    81 
    82 VMatrix NN = Group("ConcatRows", For(2,neighbour.num,VMatrix(Real j)
    83 {
    84   Mat2VMat(SubCol(neighbour::_.neighbour,[[j]]))   
    85 }));
    86 
    87 
    88 VMatrix NNSC = MatQuery::SortAndCount(NN);
    89 
    90 Matrix candidates.dens =
    91 {
    92   VMatrix aux = RPow(SubCol(NNSC,[[2]]),-1) $* Mat2VMat(candidates.isInDomain);
    93   Real sum.f = VMatSum(aux);
    94   VMatrix w = (aux / sum.f);
    95   VMat2Mat(SubCol(NNSC,[[1]]) | w )
    96 };
    97 
    98 Matrix candidates.dens.sort =
    99 {
    100   Matrix perm = Sort(candidates.dens,[[2]]);
    101   SubRow(PivotByRows(candidates.dens,perm),
    102          Range(numOutDomain+1,candidate.num,1))
    103 };
    104 
    105 Matrix candidates.dens.cum =
    106 {
    107   Matrix dens = SubCol(candidates.dens.sort,[[2]]);
    108   Matrix cum = DifEq(1/(1-B),dens);
    109   candidates.dens.sort | cum
    110 };
    111 
    112 Set workSpace =  { AlgLib.Interp.Vector.CreateWorkSpace(
    113   Text interpMet = "LinearSpline",
    114   Matrix x = Tra(SubCol(candidates.dens.cum,[[3]])),
    115   Matrix F = Tra(DifEq(1/(1-B),Constant(Columns(x),1,1)))) };
    116 
    117 Code fun = AlgLib.Interp.Vector.LoadWorkSpace(workSpace);
    118 
    119 Set candidates.extract = For(1,size,Real(Real k)
    120 {
    121   Real pos =
    122     Min(candidate.num-numOutDomain,
    123     Max(1,Round(MatDat(fun(0,Rand(0,1)),1,1))));
    124   MatDat(candidates.dens.sort,pos,1)
     25    Real total_iter := total_iter+1;
     26    Real iter := iter+1;
     27    Real lambda = Rand(-3*deg,3*deg);
     28    Polyn p := a+U*lambda;
     29    ok := IsStationary(p)
     30  });
     31  Matrix If(ok,
     32    sample := sample << SubCol(Tra(PolMat(p,deg+1,1)),rng_col));
     33  Real total_success := total_success +ok
    12534});
    12635
    127 Matrix sample = SubRow(candidates, candidates.extract);
    128 
    129 Set sample.showFaces = getFaces(sample);
    130 
    131 Matrix sample.isInDomain = SetCol(For(1,size,Real(Real i)
     36Set showFaces = SetConcat(For(1,deg-1,Set(Real i)
    13237{
    133   isInDomain(SubRow(sample,[[i]]))
     38  For(i+1,deg,Matrix(Real j)
     39  {
     40    Matrix face = SubCol(sample,[[i,j]]);
     41    Eval("Face_"<<i+"_"<<j+"=face")
     42  })
    13443}));
    135 
    136 Real sample.isInDomain.ratio = MatAvr(sample.isInDomain);
    137 
    138 /* */

  • tol_tests/tol/ARIMA/mcmc/MTryMetropolis.tol

    r53196af rdcc06082  
    33//////////////////////////////////////////////////////////////////////////////
    44{
     5  Real useOnlySingleTry = False;
    56  Real _.n = 0;
    67  Real _.k = 0;
    78  VMatrix _.x = Rand(0,0,0,0);
    8   Real _.pi_logDens_x = ?;
    99  Real _.r = ?;
    1010  Real _.s = 1;
    … …  
    1818  Real _.num_draw = 0;
    1919  Real _.num_accept = 0;
    20   Real _.num_Q_draw_pre = 0;
    21   Real _.num_Q_draw = 0;
    2220  Real _.num_logWeight = 0;
    2321
    24   ////////////////////////////////////////////////////////////////////////////
    25   Real pi_matchesDomain(VMatrix x);
    26   ////////////////////////////////////////////////////////////////////////////
     22  VMatrix _.nu = Rand(0,0,0,0);
     23  VMatrix _.L = Rand(0,0,0,0);  // Cov = L*LŽ
     24  VMatrix _.Li = Rand(0,0,0,0); // L*Li = I
     25
     26  ////////////////////////////////////////////////////////////////////////////
     27  VMatrix Q_draw(VMatrix x)
     28  ////////////////////////////////////////////////////////////////////////////
     29  {
     30    VMatrix e = Gaussian(_.n,1,0,1);
     31     _.nu + _.L*e*_.s
     32  };
     33
     34  ////////////////////////////////////////////////////////////////////////////
     35  Real Q_logDens(VMatrix x, VMatrix y)
     36  ////////////////////////////////////////////////////////////////////////////
     37  {
     38    VMatrix e = _.Li*(y-x);
     39    -0.5*VMatDat(MtMSqr(e),1,1)
     40  };
     41
     42  ////////////////////////////////////////////////////////////////////////////
     43  Real pi_logDens_aprox(VMatrix x)
     44  ////////////////////////////////////////////////////////////////////////////
     45  {
     46    VMatrix e = _.Li*(x-_.nu);
     47    -0.5*VMatDat(MtMSqr(e),1,1)   
     48  };
     49
     50  ////////////////////////////////////////////////////////////////////////////
     51  Real pi_matchesDomain(VMatrix x)
     52  ////////////////////////////////////////////////////////////////////////////
     53  {
     54    True
     55  };
    2756
    2857  ////////////////////////////////////////////////////////////////////////////
    … …  
    3160
    3261  ////////////////////////////////////////////////////////////////////////////
    33   Real Q_logDens(VMatrix x, VMatrix y) { 1 };
    34   ////////////////////////////////////////////////////////////////////////////
    35 
    36   ////////////////////////////////////////////////////////////////////////////
    37   VMatrix Q_draw_pre(VMatrix x);
    38   ////////////////////////////////////////////////////////////////////////////
    39 
    40   ////////////////////////////////////////////////////////////////////////////
    41   VMatrix Q_draw(VMatrix x)
    42   ////////////////////////////////////////////////////////////////////////////
    43   {
    44   //WriteLn("TRACE Q_draw 1 ");
    45     Real _.num_Q_draw := _.num_Q_draw + 1;
    46     Q_draw_pre(x)
    47 /*
    48     VMatrix x_ = If(!_.constrainedDomain,
    49     {
    50     //WriteLn("TRACE Q_draw 2.1 ");
    51       Real _.num_Q_draw_pre := _.num_Q_draw_pre + 1;
    52       Q_draw_pre(x)
    53     },
    54     {
    55     //WriteLn("TRACE Q_draw 2.2 ");
    56       VMatrix y = Rand(0,0,0,0);
    57       Real ok = False;
    58       Real iter = 1;
    59       While(And(!ok, iter<=_.maxDomainSearchIter),
    60       {
    61         VMatrix y := Q_draw_pre(x);
    62         Real ok := pi_matchesDomain(y);
    63         Real _.num_Q_draw_pre := _.num_Q_draw_pre + 1;
    64         Real iter := iter +1
    65       });
    66     //WriteLn("TRACE Q_draw 2.3 iter="<<iter+ " ok="<<ok);
    67       If(ok,y,
    68       {
    69         WriteLn("Cannot find a valid draw after "<<iter+" iterations","W");
    70         x
    71       })
    72     });
    73   //WriteLn("TRACE Q_draw 3 x_=\n"<<Matrix Tra(VMat2Mat(x_)));
    74     x_
    75 */
    76   };
    77 
    78   ////////////////////////////////////////////////////////////////////////////
    79   Real logWeight(VMatrix x, VMatrix y);
    80   ////////////////////////////////////////////////////////////////////////////
     62  Real logLambda(VMatrix x, VMatrix y)
     63  ////////////////////////////////////////////////////////////////////////////
     64  {
     65    -Q_logDens(x,y)
     66  };
     67
     68  ////////////////////////////////////////////////////////////////////////////
     69  Real logWeight(VMatrix x, VMatrix y)
     70  ////////////////////////////////////////////////////////////////////////////
     71  {
     72    pi_logDens(x)
     73  }
    8174
    8275  ////////////////////////////////////////////////////////////////////////////
    … …  
    9689    });
    9790    Real _.acceptRatio := _.num_accept / _.num_draw;
    98     Real Case(
     91    Real If(useOnlySingleTry, Case(
     92    _.k>1, { False },
     93    _.acceptRatio<_.acceptRatioTarget/1.01,
     94    {
     95      Real _.s := _.s/_.sFactor
     96    },
     97    _.acceptRatio>_.acceptRatioTarget*1.01,
     98    {
     99      Real _.s := _.s*_.sFactor
     100    }),
     101    Case(
    99102    _.acceptRatio<0.50,
    100103    {
    … …  
    106109      If(_.k>1,     Real _.k := _.k-1);
    107110      If(_.s<1.E+3, Real _.s := _.s*_.sFactor)
    108     });
    109 
    110 /*
    111     Real Case(
    112     _.k>1, { False },
    113     _.acceptRatio<_.acceptRatioTarget/1.01,
    114     {
    115       Real _.s := _.s/_.sFactor
    116     },
    117     _.acceptRatio>_.acceptRatioTarget*1.01,
    118     {
    119       Real _.s := _.s*_.sFactor
    120     });
    121 */
    122   //Real If(_.accept, { Real doAfterDraw(0) });
     111    }));
    123112    Real doAfterDraw(0);
    124113    _.accept
    125114  };
    126 
    127115
    128116  ////////////////////////////////////////////////////////////////////////////
    … …  
    133121    VMatrix y = Q_draw(_.x);
    134122    Real pi_logDens_y = pi_logDens(y);
    135     Real _.pi_logDens_x := pi_logDens(_.x);
     123    Real pi_logDens_x = pi_logDens(_.x);
    136124    Real _.r := Min(1, If(_.simetricQ,
    137125    {
    138       Exp(pi_logDens_y-_.pi_logDens_x)
     126      Exp(pi_logDens_y-pi_logDens_x)
    139127    },
    140128    {
    141       Exp((  pi_logDens_y + Q_logDens(  y, _.x)) -
    142           (_.pi_logDens_x + Q_logDens(_.x,   y)) )
     129      Exp((pi_logDens_y + Q_logDens(  y, _.x)) -
     130          (pi_logDens_x + Q_logDens(_.x,   y)) )
    143131    }));
    144132    Real accept_reject(y);
    … …  
    146134  };
    147135
     136
     137  ////////////////////////////////////////////////////////////////////////////
     138  VMatrix draw_multiple_try(Real void)
     139  ////////////////////////////////////////////////////////////////////////////
     140  {
     141    If(_.k==1, draw_single_try(void), {
    148142    VMatrix lwYjx = Rand(0,0,0,0);
    149143    VMatrix gYjx = Rand(0,0,0,0);
    … …  
    153147    Real lwXjy_max = ?;
    154148    Real gXjy_sum = ?;
    155 
    156   ////////////////////////////////////////////////////////////////////////////
    157   VMatrix draw_multiple_try(Real void)
    158   ////////////////////////////////////////////////////////////////////////////
    159   {
    160     If(_.k==1, draw_single_try(void), {
    161149    Real _.num_draw := _.num_draw + 1;
    162150    VMatrix J = DifEq(1/(1-B),Constant(_.k,1,1));
    … …  
    225213    _.x
    226214  }) }
     215
    227216};
    228217
     218
     219
    229220//////////////////////////////////////////////////////////////////////////////
    230 Class @MTryMetropolisGeneric : @MTryMetropolis
    231 //////////////////////////////////////////////////////////////////////////////
    232 {
    233 
    234   ////////////////////////////////////////////////////////////////////////////
    235   Real logLambda(VMatrix x, VMatrix y);
    236   ////////////////////////////////////////////////////////////////////////////
    237 
    238   ////////////////////////////////////////////////////////////////////////////
    239   Real logWeight(VMatrix x, VMatrix y)
    240   ////////////////////////////////////////////////////////////////////////////
    241   {
    242     pi_logDens(x)+Q_logDens(x,y)+logLambda(x,y)
    243   }
    244 };
    245 
    246 //////////////////////////////////////////////////////////////////////////////
    247 Class @MTryMetropolisSymmetric : @MTryMetropolis
    248 //////////////////////////////////////////////////////////////////////////////
    249 {
    250   Real _.simetricQ = True;
    251   Real _.constrainedDomain = False;
    252 
    253   ////////////////////////////////////////////////////////////////////////////
    254   Real logWeight(VMatrix x, VMatrix y)
    255   ////////////////////////////////////////////////////////////////////////////
    256   {
    257     pi_logDens(x)
    258   }
    259 };
    260 
    261 //////////////////////////////////////////////////////////////////////////////
    262 Class @MTryMetropolisSymmetricNormal : @MTryMetropolisSymmetric
    263 //////////////////////////////////////////////////////////////////////////////
    264 {
    265   VMatrix _.L;  // Cov = L*LŽ
    266 
    267   ////////////////////////////////////////////////////////////////////////////
    268   Real pi_matchesDomain(VMatrix x) { True };
    269   ////////////////////////////////////////////////////////////////////////////
    270 
    271   ////////////////////////////////////////////////////////////////////////////
    272   VMatrix Q_draw_pre(VMatrix x)
    273   ////////////////////////////////////////////////////////////////////////////
    274   {
    275     VMatrix e = Gaussian(_.n,1,0,1);
    276     x + _.L*e*_.s
    277   }
    278 };
    279 
    280 //////////////////////////////////////////////////////////////////////////////
    281 Class @MTryMetropolisTruncNormal : @MTryMetropolis
    282 //////////////////////////////////////////////////////////////////////////////
    283 {
    284   Real _.simetricQ = False;
    285   Real _.constrainedDomain = True;
    286   VMatrix _.nu = Rand(0,0,0,0);
    287   VMatrix _.L = Rand(0,0,0,0);  // Cov = L*LŽ
    288   VMatrix _.Li = Rand(0,0,0,0); // L*Li = I
    289   Real _.do_cum = False;
    290   Real _.pi_logScale = ?;
    291   Real _.iter_cum = 0;
    292   Real _.pi_cum = 0;
    293   VMatrix _.nu_cum = Rand(0,0,0,0);
    294   VMatrix _.var_cum = Rand(0,0,0,0);
    295 
    296   ////////////////////////////////////////////////////////////////////////////
    297   VMatrix Q_draw_pre(VMatrix x)
    298   ////////////////////////////////////////////////////////////////////////////
    299   {
    300     VMatrix e = Gaussian(_.n,1,0,1);
    301      _.nu + _.L*e*_.s
    302   };
    303   ////////////////////////////////////////////////////////////////////////////
    304   Real Q_logDens(VMatrix x, VMatrix y)
    305   ////////////////////////////////////////////////////////////////////////////
    306   {
    307     VMatrix e = _.Li*(y-_.nu);
    308     -0.5*VMatDat(MtMSqr(e),1,1)
    309   };
    310 
    311   ////////////////////////////////////////////////////////////////////////////
    312   Real cumulate(VMatrix x, Real x_pi_logDens, Real forze)
    313   ////////////////////////////////////////////////////////////////////////////
    314   {
    315     Real Case(
    316     !_.do_cum,
    317     {
    318     //WriteLn("TRACE @MTryMetropolisTruncNormal::cumulate 2.1");
    319       False 
    320     },
    321     VMatMax(Abs(_.x-_.nu))==0,
    322     {
    323     //WriteLn("TRACE @MTryMetropolisTruncNormal::cumulate 2.2");
    324       False 
    325     },
    326     x_pi_logDens <= -TheMaxAbsValue,
    327     {
    328     //WriteLn("TRACE @MTryMetropolisTruncNormal::cumulate 2.3");
    329       False 
    330     },
    331     IsUnknown(_.pi_logScale),
    332     {
    333     //WriteLn("TRACE @MTryMetropolisTruncNormal::cumulate 2.4");
    334       Real _.iter_cum := _.iter_cum +1;
    335       Real _.pi_logScale := x_pi_logDens;
    336       Real _.pi_cum := 1;
    337       VMatrix _.nu_cum := x;
    338       VMatrix _.var_cum := (x-_.nu)^2;
    339     //WriteLn("TRACE @MTryMetropolisTruncNormal::cumulate 2.4 _.var_cum=\n"<<Tra(VMat2Mat(_.var_cum)));
    340       True
    341     },
    342     x_pi_logDens < _.pi_logScale,
    343     {
    344     //WriteLn("TRACE @MTryMetropolisTruncNormal::cumulate 2.5");
    345       Real _.iter_cum := _.iter_cum +1;
    346       Real x_pi = Exp(x_pi_logDens-_.pi_logScale);
    347       Real _.pi_cum := _.pi_cum + x_pi;
    348       VMatrix _.nu_cum := _.nu_cum + x*x_pi;
    349       VMatrix _.var_cum := _.var_cum + ((x-_.nu)^2)*x_pi;
    350     //WriteLn("TRACE @MTryMetropolisTruncNormal::cumulate 2.5 _.var_cum=\n"<<Tra(VMat2Mat(_.var_cum)));
    351       True
    352     },
    353     x_pi_logDens > _.pi_logScale,
    354     {
    355     //WriteLn("TRACE @MTryMetropolisTruncNormal::cumulate 2.6");
    356       Real _.iter_cum := _.iter_cum +1;
    357       Real coef = Exp(_.pi_logScale-x_pi_logDens);
    358       Real _.pi_logScale := x_pi_logDens;
    359       Real _.pi_cum := _.pi_cum*coef + 1;
    360       VMatrix _.nu_cum := _.nu_cum*coef + x;
    361       VMatrix _.var_cum := _.var_cum*coef + ((x-_.nu)^2);
    362     //WriteLn("TRACE @MTryMetropolisTruncNormal::cumulate 2.6 _.var_cum=\n"<<Tra(VMat2Mat(_.var_cum)));
    363       True
    364     },
    365     1==1,
    366     {
    367     //WriteLn("TRACE @MTryMetropolisTruncNormal::cumulate 2.7");
    368       False
    369     })
    370   };
    371 
    372   ////////////////////////////////////////////////////////////////////////////
    373   Real logWeight(VMatrix x, VMatrix y)
    374   ////////////////////////////////////////////////////////////////////////////
    375   {
    376     Real y_Q_logDens = Q_logDens(x,y);
    377     Real x_Q_logDens = Q_logDens(y,x);
    378     Real x_pi_logDens = pi_logDens(x);
    379     Real cumulate(x,x_pi_logDens,False);
    380     x_pi_logDens + 0.5*(y_Q_logDens-x_Q_logDens)
    381   };
    382 
    383   ////////////////////////////////////////////////////////////////////////////
    384   Real initialize_Q(Real void)
    385   ////////////////////////////////////////////////////////////////////////////
    386   {
    387     Real _.pi_logDens_x := pi_logDens(_.x);
    388     VMatrix _.nu := _.x;
    389     VMatrix _.L  := Eye(_.n);
    390     VMatrix _.Li := Eye(_.n);
    391     True
    392   };
    393 
    394   ////////////////////////////////////////////////////////////////////////////
    395   Real initialize_cum(Real void)
    396   ////////////////////////////////////////////////////////////////////////////
    397   {
    398     Real _.do_cum := True;
    399     Real _.pi_cum := 0;
    400     VMatrix _.nu_cum := Rand(_.n,1,0,0);
    401     VMatrix _.var_cum := Rand(_.n,1,0,0);
    402     Real cumulate(_.x, pi_logDens(_.x), True);
    403     True
    404   };
    405 
    406   ////////////////////////////////////////////////////////////////////////////
    407   Real doAfterDraw_cum(Real void)
    408   ////////////////////////////////////////////////////////////////////////////
    409   {
    410     Case(
    411     Or(!_.do_cum, !_.iter_cum),
    412     {
    413       VMatrix _.nu := _.x;
    414       False 
    415     },
    416     _.pi_logDens_x > _.pi_logScale,
    417     {
    418       VMatrix _.nu := _.x;
    419       VMatrix var = _.var_cum /_.pi_cum;
    420       VMatrix _.L  := Eye(_.n,_.n,0,RPow(var,+1/2));
    421       VMatrix _.Li := Eye(_.n,_.n,0,RPow(var,-1/2));
    422       True
    423     },
    424     1==1,
    425     {
    426       VMatrix _.nu := _.nu_cum /_.pi_cum;
    427       VMatrix var = _.var_cum /_.pi_cum;
    428       VMatrix _.L  := Eye(_.n,_.n,0,RPow(var,+1/2));
    429       VMatrix _.Li := Eye(_.n,_.n,0,RPow(var,-1/2));
    430       True
    431     })
    432   }
    433 };
    434 
    435 
    436 //////////////////////////////////////////////////////////////////////////////
    437 Class @MTryMetropolis_ARIMA : @MTryMetropolisTruncNormal
     221Class @MTryMetropolis_ARIMA : @MTryMetropolis
    438222//////////////////////////////////////////////////////////////////////////////
    439223{
    … …  
    495279
    496280  ////////////////////////////////////////////////////////////////////////////
    497   Real initialize(Set arima, Real doCum)
     281  Real initialize(Set arima)
    498282  ////////////////////////////////////////////////////////////////////////////
    499283  {
    … …  
    509293    Real _.acceptRatioTarget := 0.234  + (0.5-0.234)*Exp(1-_.n);
    510294    Real _.k := Ceil(Sqrt(_.n));
    511     Real initialize_Q(0);
    512     Real If(doCum, initialize_cum(0));
    513295    ok
    514296  };
    … …  
    519301     VMatrix z,
    520302     Real sigma,
    521      Set arima,
    522      Real doCum)
     303     Set arima)
    523304  ////////////////////////////////////////////////////////////////////////////
    524305  {
    … …  
    528309      Real _.sigma = sigma
    529310    ]];
    530     Real new::initialize(arima, doCum);
     311    Real new::initialize(arima);
    531312    new
    532313  };
    … …  
    538319    Set If(_.accept, _.arima := col2arima(_.x));
    539320    Real setArima(void)
    540 //  Real doAfterDraw_cum(void)
    541321  };
    542322
    … …  
    568348  };
    569349
    570   Real _.simetricQ = True;
    571   ////////////////////////////////////////////////////////////////////////////
    572   VMatrix Q_draw(VMatrix x)
    573   ////////////////////////////////////////////////////////////////////////////
    574   {
    575     VMatrix e = Gaussian(_.n,1,0,1);
    576     y = x + _.L*e*_.s;
    577   //WriteLn("TRACE Q_draw y=\n"<<Matrix Tra(VMat2Mat(y)));
    578     y
    579   };
    580   ////////////////////////////////////////////////////////////////////////////
    581   Real logWeight(VMatrix x, VMatrix y)
    582   ////////////////////////////////////////////////////////////////////////////
    583   {
    584     pi_logDens(x)
    585   };
    586350
    587351  ////////////////////////////////////////////////////////////////////////////

  • tol_tests/tol/ARIMA/mcmc/mtrymh.mcmc.tol

    r53196af rdcc06082  
    1212Real If(!ObjectExist("Matrix","modEst.0::Correlations::COV"),
    1313{
     14  VMatrix mtm_arima::_.nu := mtm_arima::_.x;
     15  VMatrix mtm_arima::_.L := Eye(mtm_arima::_.n);
     16  VMatrix mtm_arima::_.Li :=  Eye(mtm_arima::_.n);
    1417  Real mtm_arima::_.s := 0.001
    1518},
    … …  
    2326  };
    2427  Matrix Li = LTInverse(L);
     28  VMatrix mtm_arima::_.nu := mtm_arima::_.x;
    2529  VMatrix mtm_arima::_.L := Mat2VMat(L);
    2630  VMatrix mtm_arima::_.Li := Mat2VMat(Li);
Note: See TracChangeset for help on using the changeset viewer.