Changeset 4304

Show
Ignore:
Timestamp:
03/06/12 17:55:51 (6 years ago)
Author:
vdebuen
Message:

Refs #1136

Location:
tolp/OfficialTolArchiveNetwork/ArimaTools
Files:
5 modified

Legend:

Unmodified
Added
Removed
  • tolp/OfficialTolArchiveNetwork/ArimaTools/ArimaTools.tol

    r4277 r4304  
    88////////////////////////////////////////////////////////////////////////////// 
    99[[ 
     10#Require NonLinGloOpt; 
    1011//read only autodoc 
    1112Text _.autodoc.name = "ArimaTools"; 
     
    2627//private implementation 
    2728#Embed "RandArima.tol"; 
    28 #Embed "Identify.tol"; 
     29//#Embed "Identify.tol"; 
    2930#Embed "ARMA.Sampler.tol"; 
    3031//#Embed "general.tol"; 
     
    5253 
    5354 
     55 
  • tolp/OfficialTolArchiveNetwork/ArimaTools/ARMA.Sampler.tol

    r4277 r4304  
    11 
    22//////////////////////////////////////////////////////////////////////////////////////// 
    3 Class @ARMA.Sampler 
     3Class @Sampler 
    44//////////////////////////////////////////////////////////////////////////////////////// 
    55{ 
    6   @Matrix _.noise; 
    7   Text _.label; 
    8   Set _.arima; 
    9   Set _arima.y = Copy(Empty); 
     6  Set _.varName = Copy(Empty); 
     7  Real _.varNum = ?; 
     8 
    109  Set _.eval = Copy(Empty); 
    1110  Set _eval.y = Copy(Empty); 
    12   Set _.period = Copy(Empty); 
    13   Real _.varNum = ?; 
    14   Real _.factorNum = ?; 
    15   Set _.p = Copy(Empty); 
    16   Set _.q = Copy(Empty); 
    17   Set _.d = Copy(Empty); 
    18   Set _.uncVecStr = Copy(Empty); 
     11 
    1912  Real _.llh.x = -1/0; 
    2013  Real _.llh.max = -1/0; 
     
    2215  Matrix _.mcmc = Constant(0,0,?); 
    2316  Matrix _.mcmc.llh = Constant(0,0,?); 
     17  Matrix _.mcmc.prob = Constant(0,0,?); 
     18  Real _.accept.num = ?; 
     19  Real _.accept.ratio = ?; 
     20 
    2421  Real stepSize = 0.01; 
    25   Real accept.num = ?; 
    26   Real accept.ratio = ?; 
    2722  Real tunning.acceptTarget = 0.25; 
    28   Real tunning.sMax = 1E-1; 
    29   Real tunning.sMin = 1E-6; 
     23  Real tunning.sMax = 1; 
     24  Real tunning.sMin = 1E-8; 
    3025  Real tunning.lapse = 10; 
    3126  Real cacheSample = 100; 
    3227 
    33   ////////////////////////////////////////////////////////////////////////////////////// 
    34   Static @ARMA.Sampler CreateFromLabel(Matrix noise, Text label) 
    35   ////////////////////////////////////////////////////////////////////////////////////// 
    36   { 
    37     @ARMA.Sampler new =  
    38     [[ 
    39       @Matrix _.noise = [[noise]]; 
    40       Text _.label = label; 
    41       Set _.arima = GetArimaFromLabel(label) 
    42     ]]; 
    43     Real new::initialize(?); 
    44     new 
    45   }; 
    46  
    47   ////////////////////////////////////////////////////////////////////////////////////// 
    48   Static @ARMA.Sampler CreateFromArima(Matrix noise, Set arima) 
    49   ////////////////////////////////////////////////////////////////////////////////////// 
    50   { 
    51     @ARMA.Sampler new =  
    52     [[ 
    53       @Matrix _.noise = [[noise]]; 
    54       Set _.arima = DeepCopy(arima); 
    55       Text _.label = GetLabelFromArima(_.arima) 
    56     ]]; 
    57     Real new::initialize(?); 
    58     new 
    59   }; 
    60  
    61   ////////////////////////////////////////////////////////////////////////////////////// 
    62   Real initialize(Real void) 
    63   ////////////////////////////////////////////////////////////////////////////////////// 
    64   { 
    65     Set _.period := EvalSet(_.arima,Real(@ARIMAStruct f) { f->Periodicity }); 
    66     Real _.factorNum := Card(_.period); 
    67     Set _.uncVecStr := arima2vectorStruct(_.arima); 
    68     Matrix _.x := arima2vector(_.arima); 
    69     Real _.varNum := Rows(_.x); 
     28  Set evalLogLikelihood(Matrix z); 
     29  Real setEvalMethod(Real void); 
     30  Matrix get.initialPoint(Real void); 
     31 
     32  ////////////////////////////////////////////////////////////////////////////////////// 
     33  Real set.initialPoint(Matrix x) 
     34  ////////////////////////////////////////////////////////////////////////////////////// 
     35  { 
     36    Matrix _.x := x; 
    7037    Real _.llh.x := pi.logDens(_.x); 
     38    Set _.eval := _eval.y; 
    7139    Real _.llh.max := _.llh.x; 
    7240    Matrix _.mcmc := Tra(_.x); 
    7341    Matrix _.mcmc.llh := Col(_.llh.x); 
    74     Real accept.num := 1; 
    75     Real accept.ratio := 1; 
    76     True 
    77   }; 
    78   
    79   ////////////////////////////////////////////////////////////////////////////////////// 
    80   Set arima2vectorStruct(Set arima) 
    81   ////////////////////////////////////////////////////////////////////////////////////// 
    82   { 
    83     EvalSet(arima, Set(@ARIMAStruct f) 
    84     {[[ 
    85       Matrix ar = SttPol2UncVec(f->AR, f->Periodicity); 
    86       Matrix ma = SttPol2UncVec(f->MA, f->Periodicity) 
    87     ]]}) 
    88   }; 
    89  
    90   ////////////////////////////////////////////////////////////////////////////////////// 
    91   Matrix arima2vector(Set arima) 
    92   ////////////////////////////////////////////////////////////////////////////////////// 
    93   { 
    94     Tra(GetNumeric(arima2vectorStruct(arima))) 
    95   }; 
    96  
    97   ////////////////////////////////////////////////////////////////////////////////////// 
    98   Set vector2arima(Matrix vector) 
    99   ////////////////////////////////////////////////////////////////////////////////////// 
    100   { 
    101     Set uncVecStr = DeepCopy(_.uncVecStr, vector); 
    102     For(1,_.factorNum, @ARIMAStruct(Real k) 
    103     { 
    104       @Set f = [[ _.arima[k] ]]; 
    105       @Set v = [[ uncVecStr[k] ]]; 
    106       @ARIMAStruct(  
    107         $f->Periodicity, 
    108         Polyn AR = UncVec2SttPol($v[1], $f->Periodicity); 
    109         Polyn MA = UncVec2SttPol($v[2], $f->Periodicity); 
    110         $f->DIF ) 
    111     }) 
    112   }; 
    113  
    114   ////////////////////////////////////////////////////////////////////////////////////// 
    115   Matrix Q.draw(Real void) 
     42    Real _.accept.num := 1; 
     43    Real _.accept.ratio := 1; 
     44    True 
     45  }; 
     46 
     47  ////////////////////////////////////////////////////////////////////////////////////// 
     48  Matrix Q.draw(Matrix x) 
    11649  ////////////////////////////////////////////////////////////////////////////////////// 
    11750  { 
    11851    Matrix e = Gaussian(_.varNum,1,0,stepSize); 
    119     _.x + e 
    120   }; 
    121  
    122   ////////////////////////////////////////////////////////////////////////////////////// 
    123   Set evalLogLikelihood(Matrix z) 
    124   ////////////////////////////////////////////////////////////////////////////////////// 
    125   { 
    126     Set _arima.y := vector2arima(z); 
    127     ArimaTools::@RandArimaSeries::evalLogLikelihood(_arima.y, $_.noise) 
     52    x + e 
     53  }; 
     54 
     55  ////////////////////////////////////////////////////////////////////////////////////// 
     56  Real Q.logDens(Matrix x, Matrix y) 
     57  ////////////////////////////////////////////////////////////////////////////////////// 
     58  { 
     59    0 
    12860  }; 
    12961 
     
    13668  }; 
    13769 
     70  ////////////////////////////////////////////////////////////////////////////////////// 
     71  Real log_lk_target(Matrix y, Matrix grad) 
     72  ////////////////////////////////////////////////////////////////////////////////////// 
     73  { 
     74    pi.logDens(y) 
     75  }; 
    13876 
    13977  ////////////////////////////////////////////////////////////////////////////////////// 
     
    14381    Real gamma_1 = 1;// 1/(_.num_draw^0.8); 
    14482    Real gamma_2 = gamma_1; 
    145     Real factor = Exp(0.5*gamma_2*(accept.ratio-tunning.acceptTarget)); 
     83    Real factor = Exp(0.5*gamma_2*(_.accept.ratio-tunning.acceptTarget)); 
    14684    Real s = stepSize*factor; 
    14785    Real stepSize := Min(tunning.sMax,Max(tunning.sMin,s)); 
     86  //WriteLn("TRACE tunning stepSize="<<stepSize); 
    14887    True 
    14988  }; 
     
    15392  ////////////////////////////////////////////////////////////////////////////////////// 
    15493  { 
    155     Matrix y = Q.draw(?); 
     94    Matrix y = Q.draw(_.x); 
    15695    Real llh.y = pi.logDens(y); 
    157     Real a = Exp(llh.y - _.llh.x); 
     96    Real qld.xy = Q.logDens(_.x,   y); 
     97    Real qld.yx = Q.logDens(  y, _.x); 
     98    Real a = Exp(llh.y - qld.xy - _.llh.x + qld.yx); 
    15899    Real r = Rand(0,1); 
    159100    Real accept = r<a; 
     
    164105     " a="<<a+" r="<<r+" accept="<<accept); 
    165106*/ 
    166  
    167107    If(!accept, _.x, 
    168108    { 
     
    173113      Matrix _.x := y 
    174114    }); 
    175     Real accept.num := accept.num+accept; 
    176     Real accept.ratio := accept.num/numSim; 
     115    Real _.accept.num := _.accept.num+accept; 
     116    Real _.accept.ratio := _.accept.num/numSim; 
    177117    Real AppendRows(_.mcmc,Tra(_.x)); 
    178118    Real AppendRows(_.mcmc.llh,Col(_.llh.x)); 
     
    187127  { 
    188128    Real t0 = Copy(Time); 
     129    Real setEvalMethod(?); 
     130    Real set.initialPoint(get.initialPoint(?)); 
    189131    Set For(2,sampleLength,Real(Real numSim) 
    190132    { 
     
    194136        Real elapsed = Copy(Time)-t0; 
    195137        Real percent = Round(10000*numSim/sampleLength)/100; 
    196         Real accepted = Round(10000*accept.ratio)/100; 
     138        Real accepted = Round(10000*_.accept.ratio)/100; 
    197139        Real remaining = Round(100*elapsed*(sampleLength-numSim)/numSim)/100; 
    198140        WriteLn("[@ARMA.Sampler::buil_mcmc] " 
     
    205147      numSim 
    206148    }); 
     149    Matrix _.mcmc.prob := Exp(_.mcmc.llh-_.llh.max); 
    207150    True 
    208151  } 
    209  
    210152}; 
    211153 
     154//////////////////////////////////////////////////////////////////////////////////////// 
     155Class @ARMA.Sampler : @Sampler 
     156//////////////////////////////////////////////////////////////////////////////////////// 
     157{ 
     158  Static NameBlock Options =  
     159  [[ 
     160     NameBlock Evaluator = [[ 
     161       Real Levinson = 1; 
     162       Real Almagro = 2 
     163     ]]; 
     164     NameBlock Optimizer = [[ 
     165       Real NonLinGloOpt = 1; 
     166       Real Marquardt = 2 
     167     ]] 
     168  ]]; 
     169  NameBlock config = [[Real _.undefined=? ]]; 
     170 
     171  Real setDefaultConfig(Real void) 
     172  { 
     173    NameBlock config := 
     174    [[ 
     175       Real Evaluator =  
     176         ArimaTools::@ARMA.Sampler::Options::Evaluator::Levinson; 
     177       Real Optimizer =  
     178         ArimaTools::@ARMA.Sampler::Options::Optimizer::NonLinGloOpt; 
     179       Real NonLinGloOpt.Algorithm =  
     180         NonLinGloOpt::Algorithm::LN_SBPLX 
     181    ]]; 
     182    True 
     183  }; 
     184 
     185  @Matrix _.noise; 
     186  Text _.label; 
     187  Set _.arima; 
     188  Set _.degStr = Copy(Empty); 
     189  Matrix _.degVec = Constant(0,0,?); 
     190  Set _arima.y = Copy(Empty); 
     191  Set _.period = Copy(Empty); 
     192  Real _.factorNum = ?; 
     193  Real _.s = ?; 
     194  Real _.p = ?; 
     195  Real _.q = ?; 
     196  Real _.d = ?; 
     197  Real _.N = ?; 
     198  Set _.uncVecStr = Copy(Empty); 
     199  Set _.optimizer = Copy(Empty); 
     200 
     201  Set _.evaluator = Copy(Empty); 
     202  Real _.ACovDetN = ?; 
     203  Real _.WACovInvW = ?; 
     204 
     205  ////////////////////////////////////////////////////////////////////////////////////// 
     206  Static Set DegVecToDegStr(Set period, Matrix decVec) 
     207  ////////////////////////////////////////////////////////////////////////////////////// 
     208  { 
     209    For(1,Card(period), Set(Real k) 
     210    { 
     211      Real base = (k-1)*3; 
     212      @ARIMADegreeStruct 
     213      ( 
     214        period[k], 
     215        MatDat(decVec,base+1,1), 
     216        MatDat(decVec,base+2,1), 
     217        MatDat(decVec,base+3,1) 
     218      ) 
     219    }) 
     220  }; 
     221 
     222  ////////////////////////////////////////////////////////////////////////////////////// 
     223  Static Matrix DegStrToDegVec(Set degStr) 
     224  ////////////////////////////////////////////////////////////////////////////////////// 
     225  { 
     226    SetCol(SetConcat(EvalSet(degStr, Set(@ARIMADegreeStruct f) 
     227    { 
     228      [[f->ARDegree, f->MADegree, f->DIFDegree ]] 
     229    }))) 
     230  }; 
     231 
     232  ////////////////////////////////////////////////////////////////////////////////////// 
     233  Static Set ArimaToPeriod(Set arima) 
     234  ////////////////////////////////////////////////////////////////////////////////////// 
     235  { 
     236    EvalSet(arima,Real(@ARIMAStruct f) { f->Periodicity }) 
     237  }; 
     238 
     239  ////////////////////////////////////////////////////////////////////////////////////// 
     240  Static Set ArimaToDegStr(Set arima) 
     241  ////////////////////////////////////////////////////////////////////////////////////// 
     242  { 
     243    EvalSet(arima, @ARIMADegreeStruct(@ARIMAStruct f) 
     244    { 
     245      @ARIMADegreeStruct 
     246      ( 
     247        f->Periodicity, 
     248        Degree(f->AR), 
     249        Degree(f->MA), 
     250        Degree(f->DIF) 
     251      ) 
     252    }) 
     253  }; 
     254 
     255  ////////////////////////////////////////////////////////////////////////////////////// 
     256  Static Set DegStrToArima(Set degStr) 
     257  ////////////////////////////////////////////////////////////////////////////////////// 
     258  { 
     259    EvalSet(degStr, @ARIMAStruct(@ARIMADegreeStruct f) 
     260    { 
     261      @ARIMAStruct 
     262      ( 
     263        f->Periodicity, 
     264        RandStationary(f->ARDegree), 
     265        RandStationary(f->MADegree), 
     266        RandStationary(f->DIFDegree) 
     267      ) 
     268    }) 
     269  }; 
     270 
     271  ////////////////////////////////////////////////////////////////////////////////////// 
     272  Static @ARMA.Sampler CreateFromLabel(Matrix noise, Text label) 
     273  ////////////////////////////////////////////////////////////////////////////////////// 
     274  { 
     275    @ARMA.Sampler new =  
     276    [[ 
     277      @Matrix _.noise = [[noise]]; 
     278      Text _.label = Copy(label); 
     279      Set _.arima = GetArimaFromLabel(_.label); 
     280      Set _.period = @ARMA.Sampler::ArimaToPeriod(_.arima); 
     281      Set _.degStr = @ARMA.Sampler::ArimaToDegStr(_.arima); 
     282      Matrix _.decVec = @ARMA.Sampler::DegStrToDegVec(_.degStr) 
     283    ]]; 
     284    Real new::initialize(?); 
     285    new 
     286  }; 
     287 
     288  ////////////////////////////////////////////////////////////////////////////////////// 
     289  Static @ARMA.Sampler CreateFromArima(Matrix noise, Set arima) 
     290  ////////////////////////////////////////////////////////////////////////////////////// 
     291  { 
     292    @ARMA.Sampler new =  
     293    [[ 
     294      @Matrix _.noise = [[noise]]; 
     295      Set _.arima = DeepCopy(arima); 
     296      Set _.period = @ARMA.Sampler::ArimaToPeriod(_.arima); 
     297      Text _.label = GetLabelFromArima(_.arima); 
     298      Set _.degStr = @ARMA.Sampler::ArimaToDegStr(_.arima); 
     299      Matrix _.decVec = @ARMA.Sampler::DegStrToDegVec(_.degStr) 
     300       
     301    ]]; 
     302    Real new::initialize(?); 
     303    new 
     304  }; 
     305 
     306  ////////////////////////////////////////////////////////////////////////////////////// 
     307  Static @ARMA.Sampler CreateFromDegStr(Matrix noise, Set degStr) 
     308  ////////////////////////////////////////////////////////////////////////////////////// 
     309  { 
     310    @ARMA.Sampler new =  
     311    [[ 
     312      @Matrix _.noise = [[noise]]; 
     313      Set _.degStr = DeepCopy(degStr); 
     314      Set _.arima = @ARMA.Sampler::DegStrToArima(degStr); 
     315      Set _.period = @ARMA.Sampler::ArimaToPeriod(_.arima); 
     316      Text _.label = GetLabelFromArima(_.arima); 
     317      Matrix _.decVec = @ARMA.Sampler::DegStrToDegVec(_.degStr) 
     318    ]]; 
     319    Real new::initialize(?); 
     320    new 
     321  }; 
     322 
     323  ////////////////////////////////////////////////////////////////////////////////////// 
     324  Static @ARMA.Sampler CreateFromDegVec(Matrix noise, Set period, Matrix decVec) 
     325  ////////////////////////////////////////////////////////////////////////////////////// 
     326  { 
     327    @ARMA.Sampler new =  
     328    [[ 
     329      @Matrix _.noise = [[noise]]; 
     330      Set _.period = DeepCopy(period); 
     331      Matrix _.decVec = DeepCopy(decVec); 
     332      Set _.decStr = @ARMA.Sampler::DegVecToDegStr(period, degVec); 
     333      Set _.arima = @ARMA.Sampler::DegStrToArima(_.degStr); 
     334      Text _.label = GetLabelFromArima(_.arima) 
     335    ]]; 
     336    Real new::initialize(?); 
     337    new 
     338  }; 
     339 
     340  ////////////////////////////////////////////////////////////////////////////////////// 
     341  Real initialize(Real void) 
     342  ////////////////////////////////////////////////////////////////////////////////////// 
     343  { 
     344    Real setDefaultConfig(?); 
     345    Real _.N := Rows($_.noise); 
     346    Real _.factorNum := Card(_.arima); 
     347    Set _.period := EvalSet(_.arima,Real(@ARIMAStruct f) { f->Periodicity }); 
     348    Real _.s := SetMax(_.period); 
     349    Real _.p := Degree(ARIMAGetAR(_.arima)); 
     350    Real _.q := Degree(ARIMAGetMA(_.arima)); 
     351    Real _.d := Degree(ARIMAGetDIF(_.arima)); 
     352    Matrix _.degVec := SetCol(SetConcat(EvalSet(_.degStr, Set(@ARIMADegreeStruct f) 
     353    { 
     354      [[f->ARDegree, f->MADegree, f->DIFDegree ]] 
     355    }))); 
     356    Set _.varName := SetConcat(For(1,_.factorNum,Set(Real k) 
     357    { 
     358      Real s = _.period[k]; 
     359      Text prefix = "F"<<k+".P"<<s+"."; 
     360      @Set f = [[ _.arima[k] ]]; 
     361      For(1,Degree($f->AR)/s,Text(Real j) { prefix+"AR"<<(j*s) }) << 
     362      For(1,Degree($f->MA)/s,Text(Real j) { prefix+"MA"<<(j*s) }) 
     363    })); 
     364    Set _.uncVecStr := arima2vectorStruct(_.arima); 
     365    Real _.varNum := Card(_.varName); 
     366    True 
     367  }; 
     368  
     369  //////////////////////////////////////////////////////////////////////////// 
     370  Matrix get.initialPoint(Real void) 
     371  //////////////////////////////////////////////////////////////////////////// 
     372  { 
     373    arima2vector(_.arima) 
     374  }; 
     375 
     376  //////////////////////////////////////////////////////////////////////////// 
     377  Real optimizer.NonLinGloOpt.start(Real void) 
     378  //////////////////////////////////////////////////////////////////////////// 
     379  { 
     380    @NameBlock aref = [[_this]];  
     381    Set _.optimizer := @NameBlock(NonLinGloOpt::@PipeLine::undefined(?)); 
     382    Matrix $_.optimizer::x0 := $aref::_.x; 
     383    Real $_.optimizer::set_problem(NonLinGloOpt::@Problem problem = [[ 
     384      Real sign = 1; 
     385      Real neededGlobal = False; 
     386      Real id_analyticalClass = 2; 
     387      Real id_useGradient = False; 
     388      Real n = $aref::_.varNum; 
     389      Anything target = [[$aref, log_lk_target ]]; 
     390      Set inequations = Copy(Empty) 
     391    ]]); 
     392    Real $_.optimizer::add_method(NonLinGloOpt::@Method method = 
     393    [[ 
     394      Real id_algorithm = config::NonLinGloOpt.Algorithm 
     395    ]]); 
     396    Set sc = @NameBlock(($_.optimizer::methods[1])::stopCriteria); 
     397    Real $sc::absoluteTolerance_target := .001; 
     398    Real $sc::relativeTolerance_var := .0001; 
     399    Real ($_.optimizer::problems[1])::check(?); 
     400    True 
     401  }; 
     402 
     403 
     404  //////////////////////////////////////////////////////////////////////////// 
     405  Real optimizer.NonLinGloOpt.run(Real void) 
     406  //////////////////////////////////////////////////////////////////////////// 
     407  { 
     408    Real $_.optimizer::optimize(False); 
     409    Real set.initialPoint($_.optimizer::x.opt); 
     410    True     
     411  }; 
     412 
     413  //////////////////////////////////////////////////////////////////////////// 
     414  Real optimizer.Marquardt.start(Real void) 
     415  //////////////////////////////////////////////////////////////////////////// 
     416  { 
     417    Set arimaTr = Traspose(_.arima); 
     418        Set modelDef = @ModelDef(  
     419      Serie Output = MatSerSet(Tra($_.noise),C,y2000-Rows($_.noise))[1];  
     420      Real FstTransfor = 1;  
     421      Real SndTransfor = 0;  
     422      Real Period = SetMax(arimaTr[1]);  
     423      Real Constant = 0;  
     424      Polyn Dif = SetProd(arimaTr[4]);  
     425      Set AR = arimaTr[2];  
     426      Set MA = arimaTr[3];  
     427      Set Input = Copy(Empty);  
     428      Set NonLinInput = Copy(Empty)  
     429    ); 
     430    Set _.optimizer := @NameBlock(NameBlock model =  
     431    [[ 
     432      Set def = modelDef; 
     433      Set est = Copy(Empty) 
     434    ]]); 
     435    True 
     436  }; 
     437 
     438  //////////////////////////////////////////////////////////////////////////// 
     439  Real optimizer.Marquardt.run(Real void) 
     440  //////////////////////////////////////////////////////////////////////////// 
     441  { 
     442    Set $_.optimizer::est := Estimate($_.optimizer::def); 
     443    Set res = $_.optimizer::est::Definition;  
     444    Set For(1,_.factorNum,Real(Real f) 
     445    { 
     446      Polyn _.arima[f]->AR := (res->AR)[f]; 
     447      Polyn _.arima[f]->MA := (res->MA)[f]; 
     448      f 
     449    }); 
     450    Matrix v = arima2vector(_.arima); 
     451    Real set.initialPoint(v); 
     452    True 
     453  }; 
     454 
     455  //////////////////////////////////////////////////////////////////////////// 
     456  Real optimizer.start(Real void) 
     457  //////////////////////////////////////////////////////////////////////////// 
     458  { 
     459    Real setEvalMethod(?); 
     460    Real set.initialPoint(get.initialPoint(?)); 
     461    If(config::Optimizer== 
     462       ArimaTools::@ARMA.Sampler::Options::Optimizer::NonLinGloOpt, 
     463    { 
     464      optimizer.NonLinGloOpt.start(?)  
     465    }, 
     466    { 
     467      optimizer.Marquardt.start(?)  
     468    }) 
     469  }; 
     470 
     471  //////////////////////////////////////////////////////////////////////////// 
     472  Real optimizer.run(Real void) 
     473  //////////////////////////////////////////////////////////////////////////// 
     474  { 
     475    If(config::Optimizer== 
     476       ArimaTools::@ARMA.Sampler::Options::Optimizer::NonLinGloOpt, 
     477    { 
     478      optimizer.NonLinGloOpt.run(?)  
     479    }, 
     480    { 
     481      optimizer.Marquardt.run(?)  
     482    }) 
     483  }; 
     484 
     485  ////////////////////////////////////////////////////////////////////////////////////// 
     486  Set arima2vectorStruct(Set arima) 
     487  ////////////////////////////////////////////////////////////////////////////////////// 
     488  { 
     489    EvalSet(arima, Set(@ARIMAStruct f) 
     490    {[[ 
     491      Matrix ar = SttPol2UncVec(f->AR, f->Periodicity); 
     492      Matrix ma = SttPol2UncVec(f->MA, f->Periodicity) 
     493    ]]}) 
     494  }; 
     495 
     496  ////////////////////////////////////////////////////////////////////////////////////// 
     497  Matrix arima2vector(Set arima) 
     498  ////////////////////////////////////////////////////////////////////////////////////// 
     499  { 
     500    Tra(GetNumeric(arima2vectorStruct(arima))) 
     501  }; 
     502 
     503  ////////////////////////////////////////////////////////////////////////////////////// 
     504  Set vector2arima(Matrix vector) 
     505  ////////////////////////////////////////////////////////////////////////////////////// 
     506  { 
     507    Set uncVecStr = DeepCopy(_.uncVecStr, vector); 
     508    For(1,_.factorNum, @ARIMAStruct(Real k) 
     509    { 
     510      @Set f = [[ _.arima[k] ]]; 
     511      @Set v = [[ uncVecStr[k] ]]; 
     512      @ARIMAStruct(  
     513        $f->Periodicity, 
     514        Polyn AR = UncVec2SttPol($v[1], $f->Periodicity); 
     515        Polyn MA = UncVec2SttPol($v[2], $f->Periodicity); 
     516        $f->DIF ) 
     517    }) 
     518  }; 
     519 
     520  ////////////////////////////////////////////////////////////////////////////////////// 
     521  Real setEvalMethod(Real void) 
     522  ////////////////////////////////////////////////////////////////////////////////////// 
     523  { 
     524    Case( 
     525    config::Evaluator ==  
     526       ArimaTools::@ARMA.Sampler::Options::Evaluator::Levinson, 
     527    { 
     528      Set _.evaluator := @Code(ARIMALevinsonEval), 
     529      Real _.ACovDetN := 9; 
     530      Real _.WACovInvW := 10 
     531    }, 
     532    1==1, 
     533    { 
     534      Set _.evaluator := @Code(ARIMAAlmagroEval), 
     535      Real _.ACovDetN := 6; 
     536      Real _.WACovInvW := 7 
     537    }); 
     538    True 
     539  }; 
     540 
     541  ////////////////////////////////////////////////////////////////////////////////////// 
     542  Set evalLogLikelihood(Matrix z) 
     543  ////////////////////////////////////////////////////////////////////////////////////// 
     544  { 
     545    Set _arima.y := vector2arima(z); 
     546    Code eval = $_.evaluator; 
     547    Set ev = eval(_arima.y, $_.noise); 
     548    Real N_ = Rows(ev[2]); 
     549    Real ACovDetN = ev[_.ACovDetN]; 
     550    Real WACovInvW = ev[_.WACovInvW]; 
     551    Real logLikelihood = -(N_/2)*Log(2*pi*e*ACovDetN*WACovInvW/N_); 
     552    ev<<[[logLikelihood]] 
     553  }; 
     554 
     555 
     556  //////////////////////////////////////////////////////////////////////////// 
     557  Real clear(Real void) 
     558  //////////////////////////////////////////////////////////////////////////// 
     559  {  
     560    If(config::Optimizer== 
     561       ArimaTools::@ARMA.Sampler::Options::Optimizer::NonLinGloOpt, 
     562    { 
     563      If(Card(_.optimizer), $_.optimizer::clear(?)) 
     564    }); 
     565    Set _.optimizer := Copy(Empty); 
     566    True 
     567  }; 
     568 
     569  //////////////////////////////////////////////////////////////////////////// 
     570  Real __destroy(Real void) 
     571  //////////////////////////////////////////////////////////////////////////// 
     572  {  
     573    clear(void) 
     574  } 
     575 
     576 
     577 
     578}; 
     579 
  • tolp/OfficialTolArchiveNetwork/ArimaTools/Identify.tol

    r4224 r4304  
    1  
    2 //////////////////////////////////////////////////////////////////////////////////////// 
    3 NameBlock Identify = [[ 
    4 //////////////////////////////////////////////////////////////////////////////////////// 
    5  
    6 //////////////////////////////////////////////////////////////////////////////////////// 
    7 Class @Polynomial 
     1//////////////////////////////////////////////////////////////////////////////////////// 
     2Class @PolynDegreeSampler 
    83//////////////////////////////////////////////////////////////////////////////////////// 
    94{ 
    105  Real _.min; 
    116  Real _.max; 
    12   Real _.period = 1; 
    137 
    148  Set _.transition = Copy(Empty); 
    159  Set _.transition.log = Copy(Empty); 
    1610 
    17   Real lambda = 0.5; 
     11  Real lambda = 0.8; 
    1812  Real current = 0; 
    1913 
     
    3327      Set aux = For(_.min, _.max, Real(Real y) 
    3428      { 
    35         If(x==y, 0, lambda^Abs(x-y)) 
     29        lambda^Abs(x-y) 
    3630      }); 
    3731      Real sum = SetSum(aux); 
     
    6862    MatDat(_.transition.log[n+1-_.min],m+1-_.min,1) 
    6963  } 
    70    
    7164}; 
    7265 
    7366//////////////////////////////////////////////////////////////////////////////////////// 
    74 Class @Polynomial.Unitary : @Polynomial 
    75 //////////////////////////////////////////////////////////////////////////////////////// 
    76 { 
    77   Real isValid(Polyn pol) { True }; 
    78  
    79   Polyn Q.pre_draw.pol(Real degree) 
    80   { 
    81     (1-B^_.period)^degree 
    82   }; 
    83   Polyn Q.draw.pol(Real degree) 
    84   { 
    85     (1-B^_.period)^degree 
    86   } 
    87 }; 
    88  
    89 //////////////////////////////////////////////////////////////////////////////////////// 
    90 Class @Polynomial.Starionary : @Polynomial 
    91 //////////////////////////////////////////////////////////////////////////////////////// 
    92 { 
    93   Real isValid(Polyn pol) { IsStationary(pol) }; 
    94  
    95   Polyn Q.pre_draw.pol(Real degree) 
    96   { 
    97     RandStationary(degree, _.period) 
    98   }; 
    99   Polyn Q.draw.pol(Real degree) 
    100   { 
    101     RandStationary(current, _.period) 
    102 /* 
    103     Real lambda = Rand(0,1); 
    104     Q.pre_draw.pol(degree) * lambda +  
    105     lastPolyn[degree-_.min+1] * (1-lambda)  
    106 */ 
    107   } 
    108 }; 
    109  
    110 //////////////////////////////////////////////////////////////////////////////////////// 
    111 Struct @SamplerFactor 
     67Struct @FactorDegreeSampler 
    11268//////////////////////////////////////////////////////////////////////////////////////// 
    11369{ 
    11470  Real period; 
    115   @Polynomial.Starionary ar; 
    116   @Polynomial.Starionary ma; 
    117   @Polynomial.Unitary dif 
     71  @PolynDegreeSampler ar; 
     72  @PolynDegreeSampler ma; 
     73  @PolynDegreeSampler dif 
    11874}; 
    11975 
    12076//////////////////////////////////////////////////////////////////////////////////////// 
    121 Class @Sampler 
     77Class @ARIMA.Sampler : ArimaTools::@Sampler 
    12278//////////////////////////////////////////////////////////////////////////////////////// 
    12379{ 
    124   Matrix noise; 
    125   Set factor.try; //Set of @SamplerFactor 
    126   Real factorNum = ?; 
    127   Real paramNum = ?; 
    128   Set param = Copy(Empty); 
    129   Set param.name = Copy(Empty); 
    130   Set param.min = Copy(Empty); 
    131   Set param.max = Copy(Empty); 
    132  
    133   Set auxInfo = Copy(Empty); 
    134   Real auxInfo.save = True; 
    135  
    136   Real llh.last = ?; 
    137   Real llh.max = -1/0; 
    138   Real accept.num = ?; 
    139   Real accept.ratio = ?; 
    140   Real cacheSample = 100; 
    141  
    142   Set arima.n = Copy(Empty); 
    143   Set arima.m = Copy(Empty); 
    144  
    145   Set eval = Copy(Empty); 
    146   Matrix mcmc_llh = Constant(0, 0, ?); 
    147   Matrix mcmc = Constant(0, 0, ?); 
    148  
    149   Set ranking = Copy(Empty); 
    150   Set sample = Copy(Empty); 
    151   Set histogram = Copy(Empty); 
     80  @Matrix noise; 
     81  /* 
     82  Struct @ARIMADegreeRangeStruct 
     83  { 
     84    Real Periodicity, 
     85    Real MinARDegree, 
     86    Real MaxARDegree, 
     87    Real MinMADegree, 
     88    Real MaxMADegree, 
     89    Real MinDIFDegree, 
     90    Real MaxDIFDegree 
     91  }; 
     92  */ 
     93  Set _.degreeRange;  
     94  //Set of @FactorDegreeSampler 
     95  Set _.degreeSampler = Copy(Empty); 
     96  //Set of real 
     97  Set _.period = Copy(Empty); 
     98  //Set of references to each degree samplers 
     99  Set _.degRef = Copy(Empty); 
     100  //Set of @ARMA.Sampler 
     101  Set _.armaSampler = Copy(Empty); 
     102  Real _.factorNum = ?; 
     103  Set _.varMin = Copy(Empty); 
     104  Set _.varMax = Copy(Empty); 
     105 
     106  ////////////////////////////////////////////////////////////////////////////////////// 
     107  Static @ARIMA.Sampler CreateFromDegreeRange(Matrix noise, Set degreeRange) 
     108  ////////////////////////////////////////////////////////////////////////////////////// 
     109  { 
     110    @ARIMA.Sampler new =  
     111    [[ 
     112      @Matrix _.noise = [[noise]]; 
     113      Set _.degreeRange = DeepCopy(_.degreeRange) 
     114    ]]; 
     115    Real new::initialize(?); 
     116    new 
     117  }; 
    152118 
    153119  ////////////////////////////////////////////////////////////////////////////////////// 
     
    155121  ////////////////////////////////////////////////////////////////////////////////////// 
    156122  { 
    157     If(!IsUnknown(factorNum), False, { 
    158     Real factorNum := Card(factor.try); 
    159     Set For(1, factorNum, Real(Real k) 
    160     { 
    161       @Set tk = [[ factor.try[k] ]]; 
    162       @NameBlock ar = [[ $tk->ar ]]; 
    163       @NameBlock ma = [[ $tk->ma ]]; 
    164       @NameBlock dif = [[ $tk->dif ]]; 
    165       Real $ar::_.period := $tk->period;  
    166       Real $ma::_.period := $tk->period;  
    167       Real $dif::_.period := $tk->period;  
    168       Real $ar::initialice(?);  
    169       Real $ma::initialice(?);  
    170       Real $dif::initialice(?);  
    171       Text prefix = "F"<<k+".S"<<$tk->period+"."; 
    172       Set Append(param.name, [[ 
    173         prefix+"arDeg", 
    174         prefix+"maDeg", 
    175         prefix+"difDeg" ]]); 
    176       Set Append(param, [[ 
    177         $ar::current, 
    178         $ma::current, 
    179         $dif::current ]]); 
    180       Set Append(param.min, [[ 
    181         $ar::_.min, 
    182         $ma::_.min, 
    183         $dif::_.min ]]); 
    184       Set Append(param.max, [[ 
    185         $ar::_.max, 
    186         $ma::_.max, 
    187         $dif::_.max ]]); 
    188       k 
    189     }); 
    190     Set arima.n := For(1, factorNum, Set(Real k) 
    191     { 
    192       @Set tk = [[ factor.try[k] ]]; 
    193       @ARIMAStruct( 
    194         Real $tk->period;  
    195         Polyn ($tk->ar )::Q.draw.pol(($tk->ar )::_.min);  
    196         Polyn ($tk->ma )::Q.draw.pol(($tk->ma )::_.min); 
    197         Polyn ($tk->dif)::Q.draw.pol(($tk->dif)::_.min)) 
    198     }); 
    199     Real paramNum := Card(param); 
    200     Real accept.num := 0; 
     123    Real _.factorNum := Card(_.degreeRange); 
     124    Set _.period := EvalSet(_.degreeRange, Real(@ARIMADegreeRangeStruct dr) 
     125    { 
     126      dr->period 
     127    }); 
     128    Set _.degreeSampler := EvalSet(_.degreeRange, Set(@ARIMADegreeRangeStruct dr) 
     129    { 
     130      Set fds = @FactorDegreeSampler( 
     131        Real dr->period; 
     132        @PolynDegreeSampler ar = [[ 
     133          Real _.min = dr->MinARDegree, 
     134          Real _.max = dr->MaxARDegree ]]; 
     135        @PolynDegreeSampler ma = [[ 
     136          Real _.min = dr->MinMADegree, 
     137          Real _.max = dr->MaxMADegree ]]; 
     138        @PolynDegreeSampler dif = [[ 
     139          Real _.min = dr->MinDIFDegree, 
     140          Real _.max = dr->MaxDIFDegree ]] 
     141      ); 
     142      Set Append( _.degRef, [[fds->ar,fds->ma,fds->dif]]); 
     143      fds 
     144    }); 
     145    Set _.varName := SetConcat(For(1, _.factorNum, Set(Real k) 
     146    { 
     147      Text prfx = "F."<<k+"_S."<<_.degreeRange[k]->Periodicity+"_"; 
     148      SetOfText(prfx+"AR", prfx+"MA", prfx+"DIF") 
     149    })); 
     150    Set _.varMin := SetConcat(EvalSet(_.degreeRange, Set(@ARIMADegreeRangeStruct dr) 
     151    { 
     152      SetOfReal(dr->MinARDegree, dr->MinMADegree, dr->MinDIFDegree) 
     153    })); 
     154    Set _.varMax := SetConcat(EvalSet(_.degreeRange, Set(@ARIMADegreeRangeStruct dr) 
     155    { 
     156      SetOfReal(dr->MaxARDegree, dr->MaxMADegree, dr->MaxDIFDegree) 
     157    })); 
     158    Real _.varNum := Card(_.varName); 
     159    Real set.initialPoint(Constant(_.varNum,1,0)); 
    201160    True 
    202161  })}; 
    203162 
    204163  ////////////////////////////////////////////////////////////////////////////////////// 
     164  Matrix Q.draw(Matrix x) 
     165  ////////////////////////////////////////////////////////////////////////////////////// 
     166  { 
     167    For( 
     168  }; 
     169 
     170  ////////////////////////////////////////////////////////////////////////////////////// 
     171  Real Q.logDens(Matrix x, Matrix y) 
     172  ////////////////////////////////////////////////////////////////////////////////////// 
     173  { 
     174    0 
     175  }; 
     176 
     177  ////////////////////////////////////////////////////////////////////////////////////// 
    205178  Matrix draw(Real numSim) 
    206179  ////////////////////////////////////////////////////////////////////////////////////// 
     
    208181    Real f = 0; 
    209182    Real h = 0; 
    210     @NameBlock try = [[ factor.try[1][2] ]]; 
     183    @NameBlock try = [[ factor[1][2] ]]; 
    211184    While(Or(!f, $try::_.min==$try::_.max), 
    212185    { 
    213186      Real f := IntRand(1,factorNum); 
    214187      Real h := IntRand(2,4); 
    215       @NameBlock try := [[ factor.try[f][h] ]] 
     188      @NameBlock try := [[ factor[f][h] ]] 
    216189    }); 
    217190    Real j = (f-1)*3+h-1; 
     
    350323 
    351324  ////////////////////////////////////////////////////////////////////////////////////// 
    352   Static @Sampler Create(@RandArima rndArima, Matrix noise_) 
     325  Static @Sampler Create(ArimaTools::@RandArima rndArima, Matrix noise_) 
    353326  ////////////////////////////////////////////////////////////////////////////////////// 
    354327  {  
     
    356329    [[ 
    357330      Matrix noise = noise_; 
    358       Set factor.try = For(1,Card(rndArima::factor),@SamplerFactor(Real k) 
     331      Set factor = For(1,Card(rndArima::factor),@SamplerFactor(Real k) 
    359332      { 
    360333        @NameBlock f = [[(rndArima::factor)[k] ]]; 
    361334        @SamplerFactor aux = [[ 
    362335          Real period = $f::period; 
    363           @Polynomial.Starionary ar = [[ 
     336          @PolDegree ar = [[ 
    364337             Real _.min = $f::ar::min.degree(?), 
    365338             Real _.max = $f::ar::max.degree(?) ]]; 
    366           @Polynomial.Starionary ma = [[ 
     339          @PolDegree ma = [[ 
    367340             Real _.min = $f::ma::min.degree(?), 
    368341             Real _.max = $f::ma::max.degree(?) ]]; 
    369           @Polynomial.Unitary dif = [[ 
     342          @PolDegree dif = [[ 
    370343             Real _.min = $f::dif::min.degree, 
    371344             Real _.max = $f::dif::max.degree ]] ]] 
  • tolp/OfficialTolArchiveNetwork/ArimaTools/RandArima.tol

    r4277 r4304  
    44//////////////////////////////////////////////////////////////////////////////////////// 
    55{ 
    6   Real max.degree.real; 
    7   Real min.degree.real=0; 
    8   Real max.degree.complex=0; 
    9   Real min.degree.complex=0; 
    10   Real max.root.inverse = 0.90; 
    11   Real min.root.inverse = 0.25; 
    12  
    13    
     6  Real _.max.degree.real; 
     7  Real _.min.degree.real=0; 
     8  Real _.max.degree.complex=0; 
     9  Real _.min.degree.complex=0; 
     10  Real _.max.root.inverse = 0.90; 
     11  Real _.min.root.inverse = 0.25; 
     12 
    1413  ////////////////////////////////////////////////////////////////////////////////////// 
    1514  Real min.degree(Real void) 
    1615  ////////////////////////////////////////////////////////////////////////////////////// 
    1716  { 
    18     min.degree.real+2*min.degree.complex 
     17    _.min.degree.real+2*_.min.degree.complex 
    1918  }; 
    2019 
     
    2322  ////////////////////////////////////////////////////////////////////////////////////// 
    2423  { 
    25     max.degree.real+2*max.degree.complex 
     24    _.max.degree.real+2*_.max.degree.complex 
    2625  }; 
    2726 
     
    3029  ////////////////////////////////////////////////////////////////////////////////////// 
    3130  { 
    32     Real a = Rand(min.root.inverse,max.root.inverse); 
     31    Real a = Rand(_.min.root.inverse,_.max.root.inverse); 
    3332  //WriteLn("TRACE [rand.real] a="<<a); 
    3433    1-a*B^period 
     
    3938  ////////////////////////////////////////////////////////////////////////////////////// 
    4039  { 
    41     Real a = Rand(Max(0,min.root.inverse),max.root.inverse); 
     40    Real a = Rand(Max(0,_.min.root.inverse),_.max.root.inverse); 
    4241    Real b = Rand(2*a-1,1-2*a); 
    4342  //WriteLn("TRACE [rand.complex] a="<<a+" b="<<b); 
     
    4847  ////////////////////////////////////////////////////////////////////////////////////// 
    4948  { 
    50     Real degree.real = IntRand(min.degree.real, max.degree.real); 
    51     Real degree.complex = IntRand(min.degree.complex, max.degree.complex); 
     49    Real degree.real = IntRand(_.min.degree.real, _.max.degree.real); 
     50    Real degree.complex = IntRand(_.min.degree.complex, _.max.degree.complex); 
    5251  //WriteLn("TRACE [@RandStationary::rand] degree.real = "<<degree.real); 
    5352  //WriteLn("TRACE [@RandStationary::rand] degree.complex = "<<degree.complex); 
     
    7069//////////////////////////////////////////////////////////////////////////////////////// 
    7170{ 
    72   Real min.degree=0; 
    73   Real max.degree=2; 
     71  Real _.min.degree=0; 
     72  Real _.max.degree=2; 
    7473 
    7574  Polyn draw(Real period) 
    7675  { 
    77     Real deg = IntRand(min.degree, max.degree); 
     76    Real deg = IntRand(_.min.degree, _.max.degree); 
    7877    (1-B^period)^deg 
    7978  } 
     
    8584{ 
    8685  Real period; 
    87   @RandStationary ar = [[ Real max.degree.real = 3 ]]; 
    88   @RandStationary ma = [[ Real max.degree.real = 3 ]]; 
    89   @RandUnitary dif = [[ Real max.degree=2 ]] 
     86  @RandStationary ar = [[ Real _.max.degree.real = 3 ]]; 
     87  @RandStationary ma = [[ Real _.max.degree.real = 3 ]]; 
     88  @RandUnitary dif = [[ Real _.max.degree=2 ]] 
    9089}; 
    9190 
     
    9695{ 
    9796  Real sigma = 1; 
    98   Real min.data.length = 100; 
    99   Real max.data.length = 1000; 
     97  Real _.min.data.length = 100; 
     98  Real _.max.data.length = 1000; 
    10099  Set arima = Copy(Empty); //Set of @ARIMAStruct 
    101100  Polyn ar = 1; 
     
    105104  Real q = 0; 
    106105  Real d = 0; 
     106  Real N = 0; 
    107107  Matrix residuals = Constant(0,0,?); 
    108108  Matrix noise = Constant(0,0,?); 
     
    151151  { 
    152152    Real minSize = d+2*(p+q); 
    153     Real N = IntRand(Max(min.data.length,minSize),Max(minSize,max.data.length)); 
     153    Real N := IntRand(Max(_.min.data.length,minSize),Max(minSize,_.max.data.length)); 
    154154    Real M = Min(10*N, 1E6, Max(1000, 100*N)); 
    155155    Matrix a = Gaussian(M,1,0,sigma); 
     
    167167//////////////////////////////////////////////////////////////////////////////////////// 
    168168{ 
    169   Set factor = Copy(Empty); //Set of @RandArimaFactor 
     169  Set _.factor = Copy(Empty); //Set of @RandArimaFactor 
     170  Set _.degreeRange = Copy(Empty); // Set of @ARIMADegreeRangeStruct 
     171  Real _.numCases = ?; 
    170172 
    171173  ////////////////////////////////////////////////////////////////////////////////////// 
     
    182184  ////////////////////////////////////////////////////////////////////////////////////// 
    183185  { 
    184     Set factor := EvalSet(periods,@RandArimaFactor(Real period_) 
     186    Set _.factor := EvalSet(periods,@RandArimaFactor(Real period_) 
    185187    { 
    186188      @RandArimaFactor f = [[ 
    187189        Real period = period_; 
    188         @RandUnitary dif = [[ Real max.degree=If(period_==1,2,1) ]]; 
    189         @RandStationary ar = [[ Real max.degree.real = If(period_==1,4,1) ]]; 
    190         @RandStationary ma = [[ Real max.degree.real = If(period_==1,2,1) ]] 
     190        @RandUnitary dif = [[ Real _.max.degree=If(period_==1,2,1) ]]; 
     191        @RandStationary ar = [[ Real _.max.degree.real = If(period_==1,4,1) ]]; 
     192        @RandStationary ma = [[ Real _.max.degree.real = If(period_==1,2,1) ]] 
    191193      ]]; 
    192194      f 
    193195    }); 
     196    Set _.degreeRange := EvalSet(_.factor, @ARIMADegreeRangeStruct(@RandArimaFactor f) 
     197    { 
     198      @ARIMADegreeRangeStruct 
     199      ( 
     200        Real f::period, 
     201        Real f::ar::min.degree(?), 
     202        Real f::ar::max.degree(?), 
     203        Real f::ma::min.degree(?), 
     204        Real f::ma::max.degree(?), 
     205        Real f::dif::_.min.degree, 
     206        Real f::dif::_.max.degree 
     207      ) 
     208    }); 
     209    Real _.numCases := SetProd(EvalSet(_.degreeRange, Real(@ARIMADegreeRangeStruct dr) 
     210    { 
     211      (1 + dr->MaxARDegree  - dr->MinARDegree) * 
     212      (1 + dr->MaxMADegree  - dr->MinMADegree) * 
     213      (1 + dr->MaxDIFDegree - dr->MinDIFDegree) 
     214    })); 
    194215    True 
    195216  }; 
    196217 
     218 
    197219  ////////////////////////////////////////////////////////////////////////////////////// 
    198220  Set draw(Real void) 
    199221  ////////////////////////////////////////////////////////////////////////////////////// 
    200222  { 
    201     Set arima = For(1,Card(factor),Set(Real k) 
    202     { 
    203       @NameBlock f = [[ factor[k] ]]; 
     223    Set arima = For(1,Card(_.factor),Set(Real k) 
     224    { 
     225      @NameBlock f = [[ _.factor[k] ]]; 
    204226      @ARIMAStruct( 
    205227        $f::period,  
  • tolp/OfficialTolArchiveNetwork/ArimaTools/test/test_0002/check.tol

    r4277 r4304  
    22#Require MatQuery; 
    33 
    4 Real PutRandomSeed(0); 
     4Real PutRandomSeed(1788418050); 
    55Real rndSeed = GetRandomSeed(0); 
    66WriteLn("rndSeed="<<rndSeed); 
    77 
    8 NameBlock rndArima = ArimaTools::@RandArima::Create([[1,12]]); 
     8Set periodOptions = [[ 
     9  [[1]], 
     10  [[1,5]], 
     11  [[1,7]], 
     12  [[1,12]] /* */, 
     13  [[1,24]], 
     14//[[1,364]], 
     15  [[1,24,168]] /* */ 
     16]]; 
     17 
     18Real rndPeriodIdx = IntRand(1,Card(periodOptions)); 
     19Set rndPeriod = periodOptions[rndPeriodIdx]; 
     20 
     21NameBlock rndArima = ArimaTools::@RandArima::Create(rndPeriod); 
     22 
    923Set arima = rndArima::draw(?); 
     24Text label = GetLabelFromArima(arima); 
    1025 
    1126NameBlock rndArimaSer = ArimaTools::@RandArimaSeries::Create(arima,1); 
    1227 
    13  
    1428NameBlock sampler = ArimaTools::@ARMA.Sampler::CreateFromLabel( 
    1529  rndArimaSer::noise,   
    16   GetLabelFromArima(arima) 
     30  label 
    1731); 
    1832 
    19 Real sampler::build_mcmc(2000); 
     33Real sampler::config::Evaluator :=  
     34//ArimaTools::@ARMA.Sampler::Options::Evaluator::Levinson; 
     35  ArimaTools::@ARMA.Sampler::Options::Evaluator::Almagro; 
     36 
     37Real sampler::config::Optimizer :=  
     38  ArimaTools::@ARMA.Sampler::Options::Optimizer::NonLinGloOpt; 
     39//ArimaTools::@ARMA.Sampler::Options::Optimizer::Marquardt; 
     40 
     41Real t0 = Copy(Time); 
     42Real nErr0 = Copy(NError); 
     43Real nWar0 = Copy(NWarning); 
     44 
     45//Real Show(False,"ALL"); 
     46Real sampler::optimizer.start(?); 
     47Set stopCriteria = @NameBlock(($(sampler::_.optimizer)::methods[1])::stopCriteria); 
     48Real $stopCriteria::maxTime := 10; 
     49//Real $stopCriteria::maxEval := 1000; 
     50Real $stopCriteria::absoluteTolerance_target := .01; 
     51Real $stopCriteria::relativeTolerance_var := .001; 
     52Real sampler::optimizer.run(?); 
     53 
     54//Real sampler::build_mcmc(1000); 
     55 
     56Real Show(True,"ALL"); 
     57Real nErr = Copy(NError)-nErr0; 
     58Real nWar = Copy(NWarning)-nWar0; 
     59 
     60Real time = Copy(Time)-t0; 
     61 
     62Serie z.real = MatSerSet(Tra(rndArimaSer::noise),C,y2012-Rows(rndArimaSer::noise))[1]; 
     63Serie w.real = rndArimaSer::dif : z.real; 
     64Serie a.real = MatSerSet(Tra(rndArimaSer::residuals),C,y2012-Rows(rndArimaSer::residuals))[1]; 
     65Serie a.mle  = MatSerSet(Tra(sampler::_.eval::a),C,y2012-Rows(sampler::_.eval::a))[1]; 
     66 
     67Real stdev.real = StDsS(a.real); 
     68Real stdev.estim  = StDsS(a.mle); 
     69WriteLn("stdev.real="<<stdev.real); 
     70WriteLn("stdev.estim="<<stdev.estim); 
    2071 
    2172Real logLikelihood.real = rndArimaSer::eval::logLikelihood; 
     
    2475WriteLn("logLikelihood.estim="<<logLikelihood.estim); 
    2576 
    26 Serie z.real = MatSerSet(Tra(rndArimaSer::noise),C,y2012-Rows(rndArimaSer::noise))[1]; 
    27 Serie a.real = MatSerSet(Tra(rndArimaSer::residuals),C,y2012-Rows(rndArimaSer::residuals))[1]; 
    28 Serie a.mle  = MatSerSet(Tra(sampler::_.eval::a),C,y2012-Rows(sampler::_.eval::a))[1]; 
    29  
    3077Set arima.cmp = arima << sampler::_.arima; 
    3178 
    32 /* 
    33 Text label = GetLabelFromArima (arima); 
     79Real quality = Min(1, Exp(logLikelihood.estim-logLikelihood.real)); 
    3480 
    35 Set param = ARIMAToParam(arima); 
    36  
    37 Real numVar = Card(param); 
    38  
    39 Set uncVecStr = EvalSet(arima, Set(@ARIMAStruct factor) 
    40 {[[ 
    41   Matrix ar = SttPol2UncVec(factor->AR, factor->Periodicity); 
    42   Matrix ma = SttPol2UncVec(factor->MA, factor->Periodicity) 
    43 ]]}); 
    44  
    45 Matrix uncVec = Tra(GetNumeric(uncVecStr)); 
    46 Matrix uncVec_ = uncVec+Gaussian(numVar,1,0,0.1); 
    47  
    48 Set uncVecStr_ = DeepCopy(uncVecStr, uncVec_); 
    49  
    50 Set arima_ = For(1,Card(arima), @ARIMAStruct(Real k) 
    51 { 
    52   @Set factor = [[ arima[k] ]]; 
    53   @Set vector = [[ uncVecStr_[k] ]]; 
    54   @ARIMAStruct(  
    55     $factor->Periodicity, 
    56     Polyn AR = UncVec2SttPol($vector[1], $factor->Periodicity); 
    57     Polyn MA = UncVec2SttPol($vector[2], $factor->Periodicity); 
    58     $factor->DIF ) 
    59 }); 
     81WriteLn("quality="<<quality); 
    6082 
    6183 
    62 /* 
    63 NameBlock rndArimaSer = ArimaTools::@RandArimaSeries::Create(arima,1); 
    64 NameBlock sampler = ArimaTools::Identify::@Sampler::Create( 
    65   rndArima, rndArimaSer::noise); 
    66  
    67 Real sampler::build_mcmc(10000); 
    68  
    69 //Matrix drawn = sampler::draw(1); 
    70  
    71 Real logLikelihood.real = rndArimaSer::eval::logLikelihood; 
    72 Real logLikelihood.estim = sampler::llh.max; 
    73 WriteLn("logLikelihood.real="<<logLikelihood.real); 
    74 WriteLn("logLikelihood.estim="<<logLikelihood.estim); 
    75  
    76 Serie z.real = MatSerSet(Tra(rndArimaSer::noise),C,y2012-Rows(rndArimaSer::noise))[1]; 
    77 Serie w.try = MatSerSet(Tra(sampler::eval::w),C,y2012-Rows(sampler::eval::w))[1]; 
    78 Serie w.real = rndArimaSer::dif:z.real; 
    79 Serie a.real = MatSerSet(Tra(rndArimaSer::residuals),C,y2012-Rows(rndArimaSer::residuals))[1]; 
    80 Serie a.try  = MatSerSet(Tra(sampler::eval::a),C,y2012-Rows(sampler::eval::a))[1]; 
    81 //Serie a.estim = MatSerSet(Tra(rndArimaSer::eval::a),C,y2012-Rows(rndArimaSer::eval::a))[1]; 
    82  
    83 Set degree.real = { 
    84  Set aux.1 = SetConcat(EvalSet(rndArimaSer::arima,Set(@ARIMAStruct factor) 
    85  {[[ 
    86    Degree(factor->AR)/factor->Periodicity, 
    87    Degree(factor->MA)/factor->Periodicity, 
    88    Degree(factor->DIF)/factor->Periodicity 
    89  ]]})); 
    90  Set aux.2 = For(1,Card(aux.1),Real(Real k) 
    91  { 
    92    Text name = (sampler::param.name)[k]; 
    93    Eval(name+"=aux.1[k]") 
    94  }); 
    95  [[ [[ Real {frq = 1} ]]<<aux.2 ]] 
    96 }; 
    97  
    98 Set degree.cmp = degree.real << sampler::ranking; 
    9984 
    10085/* */