Changeset e613e67 in tol for tol_tests/tol


Ignore:
Timestamp:
Apr 14, 2010, 9:25:46 AM (17 years ago)
Author:
Víctor de Buen <vdebuen@…>
Branches:
1941-docker-image-to-test, devel, master
Children:
a63ee09
Parents:
901aaa9
Message:

Reestructuración del código del test para hacerlo de forma masiva

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

Location:
tol_tests/tol/ARIMA/mcmc
Files:
5 added
1 edited

Legend:

Unmodified
Added
Removed

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

    r901aaa9 re613e67  
    33// Test de pruebas del método Multiple-Try Metropolis
    44//////////////////////////////////////////////////////////////////////////////
    5 //Random seed setting
    6 Real PutRandomSeed(1514116838);
    7 //Real PutRandomSeed(0);
    8 Real rndSeed = GetRandomSeed(0);
    9 WriteLn("Current Random Seed = "<<rndSeed);
    105
    11 Set SimulArima = {
    12   Set SimulationParameters =
    13   [[
    14     Real Sigma            =   0.1; //Desviación típica de los residuos
    15     Real MinCycles        =     6;   //Mínimo de ciclos estacionales del output
    16     Real MinData          =    40;   //Mínima longitud de la serie output
    17     Real MaxData          =  3000;   //Máxima longitud de la serie output
    18     Real MaxInterruptions =     0;   //Máximo nº de interrupciones
    19     Real MaxInputs        =     0;   //Máximo nº de inputs
    20     Real MaxRegDIF        =     0;   //Máximo nº de diferencias regulares
    21     Real MaxSeaDIF        =     0;   //Máximo nº de diferencias estacionales
    22     Real MaxRegAR         =     4;   //Máximo grado AR regular
    23     Real MaxRegMA         =     4;   //Máximo grado AR estacional
    24     Real MaxSeaAR         =     2;   //Máximo grado MA regular
    25     Real MaxSeaMA         =     2;   //Máximo grado MA estacional
    26     Real RestrictLTqp     = False;
    27   //Posibles estacionalidades
    28     Set  PeriodsTable     = [[ 1,12 ]]
    29   ]];
    30   Set Include("../SimulArima/SimulArima.tol")
    31 };
    32  
    33 Serie z = SimulArima::ARIMASeries::DifNoise;
    34 Serie e = SimulArima::ARIMASeries::Residuals;
    35 VMatrix Z = Mat2VMat(Tra(SerMat(z)));
    36 Real SIGMA = SimulArima::SimulationParameters::Sigma;
    37 Set ARIMA = SimulArima::ARIMAPolynomials::ARIMA;
    38 
    39 Matrix PARAM = SetCol(ARIMAToParam(ARIMA));
    40 
    41 /* */
    42 //////////////////////////////////////////////////////////////////////////////
    43 Class @MTryMetropolis
    44 //////////////////////////////////////////////////////////////////////////////
    45 {
    46   Real _.n = 0;
    47   Real _.k = 0;
    48   VMatrix _.x = Rand(0,0,0,0);
    49   Real _.pi_logDens_x = ?;
    50   Real _.r = ?;
    51   Real _.s = 1;
    52   Real _.sFactor = 1.001;
    53   Real _.accept = ?;
    54   Real _.acceptRatio = ?;
    55   Real _.acceptRatioTarget = ?;
    56   Real _.constrainedDomain = False;
    57   Real _.simetricQ = True;
    58   Real _.maxDomainSearchIter = 100;
    59   Real _.num_draw = 0;
    60   Real _.num_accept = 0;
    61   Real _.num_Q_draw_pre = 0;
    62   Real _.num_Q_draw = 0;
    63   Real _.num_logWeight = 0;
    64 
    65   ////////////////////////////////////////////////////////////////////////////
    66   Real pi_matchesDomain(VMatrix x);
    67   ////////////////////////////////////////////////////////////////////////////
    68 
    69   ////////////////////////////////////////////////////////////////////////////
    70   Real pi_logDens(VMatrix x);
    71   ////////////////////////////////////////////////////////////////////////////
    72 
    73   ////////////////////////////////////////////////////////////////////////////
    74   Real Q_logDens(VMatrix x, VMatrix y) { 1 };
    75   ////////////////////////////////////////////////////////////////////////////
    76 
    77   ////////////////////////////////////////////////////////////////////////////
    78   VMatrix Q_draw_pre(VMatrix x);
    79   ////////////////////////////////////////////////////////////////////////////
    80 
    81   ////////////////////////////////////////////////////////////////////////////
    82   VMatrix Q_draw(VMatrix x)
    83   ////////////////////////////////////////////////////////////////////////////
    84   {
    85   //WriteLn("TRACE Q_draw 1 ");
    86     Real _.num_Q_draw := _.num_Q_draw + 1;
    87     Q_draw_pre(x)
    88 /*
    89     VMatrix x_ = If(!_.constrainedDomain,
    90     {
    91     //WriteLn("TRACE Q_draw 2.1 ");
    92       Real _.num_Q_draw_pre := _.num_Q_draw_pre + 1;
    93       Q_draw_pre(x)
    94     },
    95     {
    96     //WriteLn("TRACE Q_draw 2.2 ");
    97       VMatrix y = Rand(0,0,0,0);
    98       Real ok = False;
    99       Real iter = 1;
    100       While(And(!ok, iter<=_.maxDomainSearchIter),
    101       {
    102         VMatrix y := Q_draw_pre(x);
    103         Real ok := pi_matchesDomain(y);
    104         Real _.num_Q_draw_pre := _.num_Q_draw_pre + 1;
    105         Real iter := iter +1
    106       });
    107     //WriteLn("TRACE Q_draw 2.3 iter="<<iter+ " ok="<<ok);
    108       If(ok,y,
    109       {
    110         WriteLn("Cannot find a valid draw after "<<iter+" iterations","W");
    111         x
    112       })
    113     });
    114   //WriteLn("TRACE Q_draw 3 x_=\n"<<Matrix Tra(VMat2Mat(x_)));
    115     x_
    116 */
    117   };
    118 
    119   ////////////////////////////////////////////////////////////////////////////
    120   Real logWeight(VMatrix x, VMatrix y);
    121   ////////////////////////////////////////////////////////////////////////////
    122 
    123   ////////////////////////////////////////////////////////////////////////////
    124   Real doAfterDraw(Real void);
    125   ////////////////////////////////////////////////////////////////////////////
    126 
    127   ////////////////////////////////////////////////////////////////////////////
    128   Real accept_reject(VMatrix y)
    129   ////////////////////////////////////////////////////////////////////////////
    130   {
    131     Real _.accept := GE(_.r,Rand(0,1));
    132   //WriteLn("TRACE accept_reject _.accept=\n"<<_.accept);
    133     Real If(_.accept,
    134     {
    135       VMatrix _.x := Copy(y);
    136       Real _.num_accept := _.num_accept +1
    137     });
    138     Real _.acceptRatio := _.num_accept / _.num_draw;
    139     Real Case(
    140     _.acceptRatio<_.acceptRatioTarget/1.01,
    141     {
    142       Real _.s := _.s/_.sFactor
    143     },
    144     _.acceptRatio>_.acceptRatioTarget*1.01,
    145     {
    146       Real _.s := _.s*_.sFactor
    147     });
    148   //Real If(_.accept, { Real doAfterDraw(0) });
    149     Real doAfterDraw(0);
    150     _.accept
    151   };
     6Real numErr.0 = Copy(NError);
     7Real numWar.0 = Copy(NWarning);
     8Set Include("generate.model.tol");
     9Set Include("maxlh.tol");
     10Set Include("MTryMetropolis.tol");
    15211
    15312
    154   ////////////////////////////////////////////////////////////////////////////
    155   VMatrix draw_single_try(Real void)
    156   ////////////////////////////////////////////////////////////////////////////
    157   {
    158     Real _.num_draw := _.num_draw + 1;
    159     VMatrix y = Q_draw(_.x);
    160     Real pi_logDens_y = pi_logDens(y);
    161     Real _.pi_logDens_x := pi_logDens(_.x);
    162     Real _.r := Min(1, If(_.simetricQ,
    163     {
    164       Exp(pi_logDens_y-_.pi_logDens_x)
    165     },
    166     {
    167       Exp((  pi_logDens_y + Q_logDens(  y, _.x)) -
    168           (_.pi_logDens_x + Q_logDens(_.x,   y)) )
    169     }));
    170     Real accept_reject(y);
    171     _.x
    172   };
     13Set k.range = Range(1,Rows(PARAM),1);
    17314
    174     VMatrix lwYjx = Rand(0,0,0,0);
    175     VMatrix gYjx = Rand(0,0,0,0);
    176     Real lwYjx_max = ?;
    177     Real gYjx_sum  = ?;
    178     VMatrix lwXjy = Rand(0,0,0,0);
    179     Real lwXjy_max = ?;
    180     Real gXjy_sum = ?;
    181 
    182   ////////////////////////////////////////////////////////////////////////////
    183   VMatrix draw_multiple_try(Real void)
    184   ////////////////////////////////////////////////////////////////////////////
    185   {
    186     If(_.k==1, draw_single_try(void), {
    187     Real _.num_draw := _.num_draw + 1;
    188     VMatrix J = DifEq(1/(1-B),Constant(_.k,1,1));
    189     Set Y = For(1,_.k, VMatrix Y_Q_draw(Real j) { Q_draw(_.x) });
    190     Real calc_numerator_Yjx(Real void)
    191     {
    192       VMatrix lwYjx := Mat2VMat(SetCol(For(1,_.k, Real(Real j)
    193       {
    194         Real _.num_logWeight := _.num_logWeight+1;
    195         logWeight(Y[j],_.x)
    196       })));
    197     //WriteLn("TRACE draw lwYjx=\n"<<Matrix Tra(VMat2Mat(lwYjx)));
    198       Real lwYjx_max := VMatMax(lwYjx);
    199     //WriteLn("TRACE draw lwYjx_max=\n"<<lwYjx_max);
    200       VMatrix gYjx := Exp(lwYjx-lwYjx_max);
    201     //WriteLn("TRACE draw gYjx=\n"<<Matrix Tra(VMat2Mat(gYjx)));
    202       Real gYjx_sum := VMatSum(gYjx);
    203     //WriteLn("TRACE draw gYjx_sum=\n"<<gYjx_sum);
    204       True
    205     };
    206     Real calc_numerator_Yjx(0);
    207 
    208     VMatrix draw_y(Real void)
    209     {
    210       Matrix fYjx = VMat2Mat(gYjx / gYjx_sum);
    211     //WriteLn("TRACE draw fYjx=\n"<<Matrix Tra(fYjx));
    212       Matrix perm = Sort(VMat2Mat(lwYjx),[[1]]);
    213     //WriteLn("TRACE draw perm=\n"<<Matrix Tra(perm));
    214       Matrix fxYp = PivotByRows(fYjx,perm);
    215     //WriteLn("TRACE draw fxYp=\n"<<Matrix Tra(fxYp));
    216       Matrix FxYp = DifEq(1/(1-B),fxYp);
    217     //WriteLn("TRACE draw FxYp=\n"<<Matrix Tra(FxYp));
    218       Real u = Rand(0,1);
    219     //WriteLn("TRACE draw Constant(_.k,1,u)=\n"<<Matrix Tra(Constant(_.k,1,u)));
    220       Matrix FLTu = LT(FxYp,Constant(_.k,1,u));
    221     //WriteLn("TRACE draw FLTu=\n"<<Matrix Tra(FLTu));
    222       Real sFLTu = MatSum(FLTu);
    223       Real h = Min(_.k,1+sFLTu);
    224       Real j = MatDat(perm,h,1);
    225     //WriteLn("TRACE draw u="<<u+" sFLTu="<<sFLTu+" h="<<h+" j="<<j);
    226       Y[j]
    227     };
    228     VMatrix y = If(_.k==1,Y[1],draw_y(0));
    229   //WriteLn("TRACE draw pi_matchesDomain(y)="<<pi_matchesDomain(y));
    230     Set X = For(1,_.k-1, VMatrix X_Q_draw(Real j) { Q_draw(y) }) << [[_.x]];
    231     Real calc_denominator_Xjy(Real void)
    232     {
    233       VMatrix lwXjy := Mat2VMat(SetCol(For(1,_.k, Real(Real j)
    234       {
    235         Real _.num_logWeight := _.num_logWeight+1;
    236         logWeight(X[j],y)
    237       })));
    238     //WriteLn("TRACE draw lwXjy=\n"<<Matrix Tra(VMat2Mat(lwXjy)));
    239       Real lwXjy_max := VMatMax(lwXjy);
    240     //WriteLn("TRACE draw lwXjy_max=\n"<<lwXjy_max);
    241       VMatrix gXjy = Exp(lwXjy-lwXjy_max);
    242       Real gXjy_sum := VMatSum(gXjy);
    243     //WriteLn("TRACE draw gXjy_sum=\n"<<gXjy_sum);
    244       True
    245     };
    246     Real calc_denominator_Xjy(0);
    247 
    248     Real _.r := Min(1,Exp(lwYjx_max-lwXjy_max)*(gYjx_sum/gXjy_sum));
    249   //WriteLn("TRACE draw _.r=\n"<<_.r);
    250     Real accept_reject(y);
    251     _.x
    252   }) }
    253 };
    254 
    255 //////////////////////////////////////////////////////////////////////////////
    256 Class @MTryMetropolisGeneric : @MTryMetropolis
    257 //////////////////////////////////////////////////////////////////////////////
     15Set simulation = EvalSet(k.range, Set(Real k)
    25816{
    259 
    260   ////////////////////////////////////////////////////////////////////////////
    261   Real logLambda(VMatrix x, VMatrix y);
    262   ////////////////////////////////////////////////////////////////////////////
    263 
    264   ////////////////////////////////////////////////////////////////////////////
    265   Real logWeight(VMatrix x, VMatrix y)
    266   ////////////////////////////////////////////////////////////////////////////
    267   {
    268     pi_logDens(x)+Q_logDens(x,y)+logLambda(x,y)
    269   }
    270 };
    271 
    272 //////////////////////////////////////////////////////////////////////////////
    273 Class @MTryMetropolisSymmetric : @MTryMetropolis
    274 //////////////////////////////////////////////////////////////////////////////
    275 {
    276   Real _.simetricQ = True;
    277   Real _.constrainedDomain = False;
    278 
    279   ////////////////////////////////////////////////////////////////////////////
    280   Real logWeight(VMatrix x, VMatrix y)
    281   ////////////////////////////////////////////////////////////////////////////
    282   {
    283     pi_logDens(x)
    284   }
    285 };
    286 
    287 //////////////////////////////////////////////////////////////////////////////
    288 Class @MTryMetropolisSymmetricNormal : @MTryMetropolisSymmetric
    289 //////////////////////////////////////////////////////////////////////////////
    290 {
    291   VMatrix _.L;  // Cov = L*LŽ
    292 
    293   ////////////////////////////////////////////////////////////////////////////
    294   Real pi_matchesDomain(VMatrix x) { True };
    295   ////////////////////////////////////////////////////////////////////////////
    296 
    297   ////////////////////////////////////////////////////////////////////////////
    298   VMatrix Q_draw_pre(VMatrix x)
    299   ////////////////////////////////////////////////////////////////////////////
    300   {
    301     VMatrix e = Gaussian(_.n,1,0,1);
    302     x + _.L*e*_.s
    303   }
    304 };
    305 
    306 //////////////////////////////////////////////////////////////////////////////
    307 Class @MTryMetropolisTruncNormal : @MTryMetropolis
    308 //////////////////////////////////////////////////////////////////////////////
    309 {
    310   Real _.simetricQ = False;
    311   Real _.constrainedDomain = True;
    312   VMatrix _.nu = Rand(0,0,0,0);
    313   VMatrix _.L = Rand(0,0,0,0);  // Cov = L*LŽ
    314   VMatrix _.Li = Rand(0,0,0,0); // L*Li = I
    315   Real _.do_cum = False;
    316   Real _.pi_logScale = ?;
    317   Real _.iter_cum = 0;
    318   Real _.pi_cum = 0;
    319   VMatrix _.nu_cum = Rand(0,0,0,0);
    320   VMatrix _.var_cum = Rand(0,0,0,0);
    321 
    322   ////////////////////////////////////////////////////////////////////////////
    323   VMatrix Q_draw_pre(VMatrix x)
    324   ////////////////////////////////////////////////////////////////////////////
    325   {
    326     VMatrix e = Gaussian(_.n,1,0,1);
    327      _.nu + _.L*e*_.s
    328   };
    329   ////////////////////////////////////////////////////////////////////////////
    330   Real Q_logDens(VMatrix x, VMatrix y)
    331   ////////////////////////////////////////////////////////////////////////////
    332   {
    333     VMatrix e = _.Li*(y-_.nu);
    334     -0.5*VMatDat(MtMSqr(e),1,1)
    335   };
    336 
    337   ////////////////////////////////////////////////////////////////////////////
    338   Real cumulate(VMatrix x, Real x_pi_logDens, Real forze)
    339   ////////////////////////////////////////////////////////////////////////////
    340   {
    341     Real Case(
    342     !_.do_cum,
    343     {
    344     //WriteLn("TRACE @MTryMetropolisTruncNormal::cumulate 2.1");
    345       False 
    346     },
    347     VMatMax(Abs(_.x-_.nu))==0,
    348     {
    349     //WriteLn("TRACE @MTryMetropolisTruncNormal::cumulate 2.2");
    350       False 
    351     },
    352     x_pi_logDens <= -TheMaxAbsValue,
    353     {
    354     //WriteLn("TRACE @MTryMetropolisTruncNormal::cumulate 2.3");
    355       False 
    356     },
    357     IsUnknown(_.pi_logScale),
    358     {
    359     //WriteLn("TRACE @MTryMetropolisTruncNormal::cumulate 2.4");
    360       Real _.iter_cum := _.iter_cum +1;
    361       Real _.pi_logScale := x_pi_logDens;
    362       Real _.pi_cum := 1;
    363       VMatrix _.nu_cum := x;
    364       VMatrix _.var_cum := (x-_.nu)^2;
    365     //WriteLn("TRACE @MTryMetropolisTruncNormal::cumulate 2.4 _.var_cum=\n"<<Tra(VMat2Mat(_.var_cum)));
    366       True
    367     },
    368     x_pi_logDens < _.pi_logScale,
    369     {
    370     //WriteLn("TRACE @MTryMetropolisTruncNormal::cumulate 2.5");
    371       Real _.iter_cum := _.iter_cum +1;
    372       Real x_pi = Exp(x_pi_logDens-_.pi_logScale);
    373       Real _.pi_cum := _.pi_cum + x_pi;
    374       VMatrix _.nu_cum := _.nu_cum + x*x_pi;
    375       VMatrix _.var_cum := _.var_cum + ((x-_.nu)^2)*x_pi;
    376     //WriteLn("TRACE @MTryMetropolisTruncNormal::cumulate 2.5 _.var_cum=\n"<<Tra(VMat2Mat(_.var_cum)));
    377       True
    378     },
    379     x_pi_logDens > _.pi_logScale,
    380     {
    381     //WriteLn("TRACE @MTryMetropolisTruncNormal::cumulate 2.6");
    382       Real _.iter_cum := _.iter_cum +1;
    383       Real coef = Exp(_.pi_logScale-x_pi_logDens);
    384       Real _.pi_logScale := x_pi_logDens;
    385       Real _.pi_cum := _.pi_cum*coef + 1;
    386       VMatrix _.nu_cum := _.nu_cum*coef + x;
    387       VMatrix _.var_cum := _.var_cum*coef + ((x-_.nu)^2);
    388     //WriteLn("TRACE @MTryMetropolisTruncNormal::cumulate 2.6 _.var_cum=\n"<<Tra(VMat2Mat(_.var_cum)));
    389       True
    390     },
    391     1==1,
    392     {
    393     //WriteLn("TRACE @MTryMetropolisTruncNormal::cumulate 2.7");
    394       False
    395     })
    396   };
    397 
    398   ////////////////////////////////////////////////////////////////////////////
    399   Real logWeight(VMatrix x, VMatrix y)
    400   ////////////////////////////////////////////////////////////////////////////
    401   {
    402     Real y_Q_logDens = Q_logDens(x,y);
    403     Real x_Q_logDens = Q_logDens(y,x);
    404     Real x_pi_logDens = pi_logDens(x);
    405     Real cumulate(x,x_pi_logDens,False);
    406     x_pi_logDens + 0.5*(y_Q_logDens-x_Q_logDens)
    407   };
    408 
    409   ////////////////////////////////////////////////////////////////////////////
    410   Real initialize_Q(Real void)
    411   ////////////////////////////////////////////////////////////////////////////
    412   {
    413     Real _.pi_logDens_x := pi_logDens(_.x);
    414     VMatrix _.nu := _.x;
    415     VMatrix _.L  := Eye(_.n);
    416     VMatrix _.Li := Eye(_.n);
    417     True
    418   };
    419 
    420   ////////////////////////////////////////////////////////////////////////////
    421   Real initialize_cum(Real void)
    422   ////////////////////////////////////////////////////////////////////////////
    423   {
    424     Real _.do_cum := True;
    425     Real _.pi_cum := 0;
    426     VMatrix _.nu_cum := Rand(_.n,1,0,0);
    427     VMatrix _.var_cum := Rand(_.n,1,0,0);
    428     Real cumulate(_.x, pi_logDens(_.x), True);
    429     True
    430   };
    431 
    432   ////////////////////////////////////////////////////////////////////////////
    433   Real doAfterDraw_cum(Real void)
    434   ////////////////////////////////////////////////////////////////////////////
    435   {
    436     Case(
    437     Or(!_.do_cum, !_.iter_cum),
    438     {
    439       VMatrix _.nu := _.x;
    440       False 
    441     },
    442     _.pi_logDens_x > _.pi_logScale,
    443     {
    444       VMatrix _.nu := _.x;
    445       VMatrix var = _.var_cum /_.pi_cum;
    446       VMatrix _.L  := Eye(_.n,_.n,0,RPow(var,+1/2));
    447       VMatrix _.Li := Eye(_.n,_.n,0,RPow(var,-1/2));
    448       True
    449     },
    450     1==1,
    451     {
    452       VMatrix _.nu := _.nu_cum /_.pi_cum;
    453       VMatrix var = _.var_cum /_.pi_cum;
    454       VMatrix _.L  := Eye(_.n,_.n,0,RPow(var,+1/2));
    455       VMatrix _.Li := Eye(_.n,_.n,0,RPow(var,-1/2));
    456       True
    457     })
    458   }
    459 };
    460 
    461 
    462 //////////////////////////////////////////////////////////////////////////////
    463 Class @MTryMetropolis_ARIMA : @MTryMetropolisTruncNormal
    464 //////////////////////////////////////////////////////////////////////////////
    465 {
    466   VMatrix _.z;
    467   Real _.sigma;
    468   Set _.arima = Copy(Empty);
    469   Set _.arma = Copy(Empty);
    470   Polyn  _.ar = 1;
    471   Polyn  _.ma = 1;
    472   NameBlock _.almagro = [[ Real void=? ]];
    473   VMatrix _.u0 = Rand(0,0,0,0);
    474   VMatrix _.z0 = Rand(0,0,0,0);
    475   VMatrix _.a0 = Rand(0,0,0,0);
    476  
    477   ////////////////////////////////////////////////////////////////////////////
    478   VMatrix arima2col(Set arima)
    479   ////////////////////////////////////////////////////////////////////////////
    480   {
    481   //Mat2VMat(SetCol(ARIMAToParam(arima)))
    482     Set arma = EvalSet(arima, Set(ARIMAStruct f)
    483     {
    484       SetOfPolyn(1-f->AR, 1-f->MA)
    485     });
    486     Mat2VMat(Tra(GetNumeric(arma)))
    487   };
    488  
    489   ////////////////////////////////////////////////////////////////////////////
    490   Set col2arima(VMatrix x)
    491   ////////////////////////////////////////////////////////////////////////////
    492   {
    493 //  ParamToARIMA(_.arima, MatSet(Tra(VMat2Mat(x)))[1])
    494     Set arma = DeepCopy(_.arma, VMat2Mat(x));
    495     Set arima =  DeepCopy(_.arima);
    496     Set For(1,Card(_.arima), Real(Real k)
    497     {
    498       Polyn arima[k]->AR := 1-arma[k][1];
    499       Polyn arima[k]->MA := 1-arma[k][2];
    500       True
    501     });
    502     arima
    503   };
    504 
    505   ////////////////////////////////////////////////////////////////////////////
    506   Real setArima(Real void)
    507   ////////////////////////////////////////////////////////////////////////////
    508   {
    509     Real If(_.accept,
    510     {
    511       Polyn  _.ar := ARIMAGetAR(_.arima);
    512       Polyn  _.ma := ARIMAGetMA(_.arima);
    513       NameBlock _.almagro := ARMAProcess::Eval.Almagro(_.ar,_.ma,_.z,_.sigma);
    514       True
    515     });
    516     VMatrix _.u0 := _.almagro::Draw.U_cond_Z(0);
    517     VMatrix _.z0 := _.almagro::Get.Z0(_.u0);
    518     VMatrix _.a0 := _.almagro::Get.A0(_.u0);
    519     True
    520   };
    521 
    522   ////////////////////////////////////////////////////////////////////////////
    523   Real initialize(Set arima, Real doCum)
    524   ////////////////////////////////////////////////////////////////////////////
    525   {
    526     Set _.arima := DeepCopy(arima);
    527     Set _.arma := EvalSet(_.arima, Set(ARIMAStruct f)
    528     {
    529       SetOfPolyn(1-f->AR, 1-f->MA)
    530     });
    531     VMatrix _.x := arima2col(_.arima);
    532     Real _.accept := True;
    533     Real ok = setArima(0);
    534     Real _.n := VRows(_.x);
    535     Real _.acceptRatioTarget := 0.234  + (0.5-0.234)*Exp(1-_.n);
    536     Real _.k := Ceil(Sqrt(_.n));
    537     Real initialize_Q(0);
    538     Real If(doCum, initialize_cum(0));
    539     ok
    540   };
    541  
    542   ////////////////////////////////////////////////////////////////////////////
    543   //Inicialización de la cadena
    544   Static @MTryMetropolis_ARIMA New(
    545      VMatrix z,
    546      Real sigma,
    547      Set arima,
    548      Real doCum)
    549   ////////////////////////////////////////////////////////////////////////////
    550   {
    551     @MTryMetropolis_ARIMA new =
    552     [[
    553       VMatrix _.z = z;
    554       Real _.sigma = sigma
    555     ]];
    556     Real new::initialize(arima, doCum);
    557     new
    558   };
    559  
    560   ////////////////////////////////////////////////////////////////////////////
    561   Real doAfterDraw(Real void)
    562   ////////////////////////////////////////////////////////////////////////////
    563   {
    564     Set If(_.accept, _.arima := col2arima(_.x));
    565     Real setArima(void)
    566 //  Real doAfterDraw_cum(void)
    567   };
    568 
    569   ////////////////////////////////////////////////////////////////////////////
    570   Real pi_matchesDomain(VMatrix x)
    571   ////////////////////////////////////////////////////////////////////////////
    572   {
    573     Set arima = col2arima(x);
    574     Polyn ar = ARIMAGetAR(arima);
    575     Polyn ma = ARIMAGetMA(arima);
    576     And(IsStationary(ar), IsStationary(ma))
    577   };
    578 
    579   ////////////////////////////////////////////////////////////////////////////
    580   Real pi_logDens_(Set arima)
    581   ////////////////////////////////////////////////////////////////////////////
    582   {
    583     Polyn ar = ARIMAGetAR(arima);
    584     Polyn ma = ARIMAGetMA(arima);
    585     If(Or(!IsStationary(ar), !IsStationary(ma)), -TheMaxAbsValue,
    586       _.almagro::LogLH.Z_cond_U(ar,ma,_.z0,_.a0) )
    587   };
    588  
    589   ////////////////////////////////////////////////////////////////////////////
    590   Real pi_logDens(VMatrix x)
    591   ////////////////////////////////////////////////////////////////////////////
    592   {
    593     pi_logDens_(col2arima(x))
    594   };
    595 
    596   Real _.simetricQ = True;
    597   ////////////////////////////////////////////////////////////////////////////
    598   VMatrix Q_draw(VMatrix x)
    599   ////////////////////////////////////////////////////////////////////////////
    600   {
    601     VMatrix e = Gaussian(_.n,1,0,1);
    602     y = x + _.L*e*_.s;
    603   //WriteLn("TRACE Q_draw y=\n"<<Matrix Tra(VMat2Mat(y)));
    604     y
    605   };
    606   ////////////////////////////////////////////////////////////////////////////
    607   Real logWeight(VMatrix x, VMatrix y)
    608   ////////////////////////////////////////////////////////////////////////////
    609   {
    610     pi_logDens(x)
    611   };
    612 
    613   ////////////////////////////////////////////////////////////////////////////
    614   Matrix move_to_average(Real maxIter)
    615   ////////////////////////////////////////////////////////////////////////////
    616   {
    617 /* * /
    618     Real pi_logDens_x = _.pi_logDens_x;
    619     VMatrix x = _.x;
    620     Matrix Y = Group("ConcatRows", For(1,maxIter,Matrix(Real iter)
    621     {
    622       VMatrix y = arima2col(RandCopyARIMA(_.arima));
    623       Real pi_logDens_y = pi_logDens(y);
    624       Real pi_logDens_y_dif = pi_logDens_y - pi_logDens_x;
    625       Real If(pi_logDens_y_dif>0,
    626       {
    627         Real coef = Exp(-pi_logDens_y_dif);
    628         VMatrix x := (x*coef + y)/(coef+1);
    629         Real pi_logDens_x := pi_logDens(x)
    630       },
    631       {
    632         Real coef = Exp(pi_logDens_y_dif);
    633         VMatrix x := (x + y*coef)/(coef+1);
    634         Real pi_logDens_x := pi_logDens(x)
    635       });
    636       Copy(VMat2Mat(Tra(x)))
    637     }));
    638 /* */
    639     Matrix pi_logDens_y = Rand(1,maxIter,0,0);
    640     Matrix Y = Group("ConcatRows", For(1,maxIter,Matrix(Real iter)
    641     {
    642       VMatrix y = arima2col(RandCopyARIMA(_.arima));
    643       Real PutMatDat(pi_logDens_y, 1, iter, pi_logDens(y));
    644       Copy(VMat2Mat(Tra(y)))
    645     }));
    646     Real pi_logDens_max = MatMax(pi_logDens_y);
    647     Matrix w = Exp(pi_logDens_y - pi_logDens_max);
    648     Matrix g = w * (1/MatSum(w));
    649     VMatrix _.x := Mat2VMat(Tra(g*Y));
    650     Real _.pi_logDens_x := pi_logDens(_.x);
    651     Real doAfterDraw(0);
    652 /* */
    653     Y
    654   }
    655 
    656 };
    657 
    658 Real PutRandomSeed(0);
    659 Real rndSeed_iv = GetRandomSeed(0);
    660 WriteLn("Current Random Seed for initial values= "<<rndSeed_iv);
    661 
    662 @MTryMetropolis_ARIMA mtm_arima = @MTryMetropolis_ARIMA::New(
    663   Z,
    664   SIGMA,
    665   RandCopyARIMA(ARIMA),
    666   False);
    667 
    668 /* * /
    669 //Ponemos un valor inicial cercano al real
    670 VMatrix mtm_arima::_.x :=
    671   mtm_arima::_.x*0.1 +
    672   mtm_arima::arima2col(ARIMA)*0.9;
    673 Real mtm_arima::doAfterDraw(0);
    674 /* */
    675 
    676 Real mtm_arima::_.k := 9; //Ceil(Sqrt(mtm_arima::_.n));
    677 
    678 Matrix mcmc = Rand(0,mtm_arima::_.n,0,0);
    679 
    680 Real mtm_arima::_.s := 0.01;
    681 
    682 /*
    683 Real m2avr.numSim = 500;
    684 Real m2avr.t0 = Copy(Time);
    685 Matrix m2avr = mtm_arima::move_to_average(m2avr.numSim);
    686 Real m2avr.time = (Copy(Time) - m2avr.t0)/m2avr.numSim;
    687 WriteLn("Time by AVR estimation "<<m2avr.time);
    688 */
    689 
    690 Real numSim = 0;
    691 
    692 Real append_draw(VMatrix draw)
    693 {
    694   If(!(numSim%100), WriteLn("
    695 ##############################################################################
    696 Simulation "<<numSim+ " _.acceptRatio="<<mtm_arima::_.acceptRatio));
    697   Matrix mcmc := mcmc << VMat2Mat(Tra(draw));
    698   Real numSim := numSim+1
    699 };
    700 
    701 /* */
    702 Real mtm_arima::_.sFactor := 1.01;
    703 Real numSim := 0;
    704 Real mh.t0 = Copy(Time);
    705 While(numSim<1000, append_draw(mtm_arima::draw_single_try(0)));
    706 Real mh.time = (Copy(Time) - mh.t0)/numSim;
    707 WriteLn("Time by MH simulation "<<mh.time);
    708 
    709 /* */
    710 Real mtm_arima::_.sFactor := 1.00;
    711 Real numSim := 0;
    712 Real mtm.t0 = Copy(Time);
    713 Real TolOprProfiler.Enabled := True;
    714 While(numSim<2500, append_draw(mtm_arima::draw_multiple_try(0)));
    715 Real TolOprProfiler.Dump("draw_multiple_try.k"<<mtm_arima::_.k+".csv");
    716 Real mtm.time = (Copy(Time) - mtm.t0)/numSim;
    717 WriteLn("Time by MTM simulation "<<mtm.time);
    718 /* */
    719 
    720 Matrix cmp = mcmc | Group("ConcatColumns", For(1,Rows(PARAM),Matrix(Real j)
    721 {
    722   Constant(Rows(mcmc),1,MatDat(PARAM,j,1))
    723 }));
    724 
    725 /* */
    726 
    727 
    728 
     17  Set aux = IncludeText(ReadFile("mtrymh.mcmc.tol"));
     18  Text name = "R"<<rndSeed+
     19              "I"<<rndSeed_iv+
     20              "L"<<arima.label+
     21              "K"<<k;
     22  Eval(name+"=aux")
     23});
Note: See TracChangeset for help on using the changeset viewer.