Changeset b7a4f8b in tol


Ignore:
Timestamp:
Dec 14, 2011, 4:20:23 PM (15 years ago)
Author:
Víctor de Buen <vdebuen@…>
Branches:
1941-docker-image-to-test, devel, master
Children:
02ea263
Parents:
58910ad
Message:

Fixes #1357

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

File:
1 edited

Legend:

Unmodified
Added
Removed

  • tol/stdlib/math/stat/models/bayesian/arima/_arma_process.tol

    r58910ad rb7a4f8b  
    1010[[
    1111  Text _.autodoc.description = "Methods related with ARMA process";
     12 
     13  ////////////////////////////////////////////////////////////////////////////
     14  Text _.autodoc.member.changeZeroCoefficient =
     15  "Change zero coefficients of a polynomial by a given value";
     16  Polyn changeZeroCoefficient(Polyn pol, Real nonZeroValue)
     17  ////////////////////////////////////////////////////////////////////////////
     18  {
     19    Matrix M = PolMat(pol, Degree(pol)+1, 1);
     20    Matrix M_ = M + Not(M)*nonZeroValue;
     21    MatPol(Tra(M_))
     22  };
    1223 
    1324  ////////////////////////////////////////////////////////////////////////////
    … …  
    95106  //WriteLn("TRACE CholeskiFactorCov("+FullName(cov)+") VRows(cov)="<<VRows(cov)+" Sparsity(cov)="<<Real(VNonNullCells(cov)/VRows(cov)^2));
    96107    VMatrix c = Drop(cov,1.E-16*VRows(cov)^2);
    97     VMatrix L = CholeskiFactor(c,"X",True);
    98   //WriteLn("TRACE CholeskiFactorCov L:"<<L);
    99     VMatrix pL = Pack(L,1);
    100   //WriteLn("TRACE CholeskiFactorCov pL:"<<pL);
    101     pL
     108    VMatrix L = CholeskiFactor(c,"X",True,False);
     109    If(VRows(L)==VRows(cov),
     110    {
     111    //WriteLn("TRACE CholeskiFactorCov L:"<<L);
     112      VMatrix pL = Pack(L,1);
     113    //WriteLn("TRACE CholeskiFactorCov pL:"<<pL);
     114      pL
     115    },
     116    {
     117      Warning("[StdLib::ARMAProcess::CholeskiFactorCov] No es posible "
     118       "aplicar CholeskiFactor a una matriz virtual no definida positiva");
     119      L
     120    })
    102121  };
    103122
    … …  
    291310  //WriteLn("TRACE Eval.Almagro Begin "+Time);
    292311  //WriteLn("TRACE Eval.Almagro(p"<<p+",q="<<q+",m="<<m);
    293     If(!p & !q, [[
     312  //Polyn   _.ar    = ARMAProcess::changeZeroCoefficient(ar,1E-5);
     313  //Polyn   _.ma    = ARMAProcess::changeZeroCoefficient(ma,1E-5);
    294314    Polyn   _.ar    = ar;
    295     Polyn   _.ma    = ma;
     315    Polyn   _.ma    = ma;
     316
     317    NameBlock neutral = [[
     318    Polyn   _.ar;
     319    Polyn   _.ma;
    296320    VMatrix _.z     = z_;
    297321    VMatrix _.z0    = Rand(0,0,0,0);
    … …  
    310334    Real _0.5_m_log2_PI_s2 = -0.5*_.m*Log(2*PI*_.s2);
    311335    Real _.logLH = _0.5_m_log2_PI_s2 - 0.5*_.m* VMatMoment(_.a,2)/_.s2
    312     ]],[[
    313     Polyn   _.ar    = ar;
    314     Polyn   _.ma    = ma;
     336    ]];
     337    NameBlock nonStationary = [[
     338    ////////////////////////////////////////////////////////////////////////////
     339    Real LogLH.Z(VMatrix z)
     340    ////////////////////////////////////////////////////////////////////////////
     341    {
     342      VMatrix _.z := z;
     343      VMatrix _.pi_z := DifEq(_.ar/_.ma,_.z);
     344      VMatrix _.G_pi_z := Tra(Tra(_.pi_z)*_.G);
     345    //Real _trc.8 = { WriteLn("TRACE Eval.Almagro Creating various "+Time); 0 };
     346      VMatrix _.u := CholeskiSolve(_.L_G, _.G_pi_z,"PtLLtP");
     347      VMatrix _.z0 := Get.Z0(_.u);
     348      VMatrix _.a0 := Get.A0(_.u);
     349      VMatrix _.v := CholeskiSolve(_.L_u, _.u,"PtL");
     350      VMatrix _.a := DifEq(_.ar/_.ma, _.z, _.z0, _.a0);
     351
     352      Real aa = MatDat(VMat2Mat( Tra(_.a)*_.a), 1,1);     
     353      Real uu = MatDat(VMat2Mat( Tra(_.u)*_.covi_u*_.u),1,1);
     354      Real _.logLH := _0.5_m_log2_PI_s2 - 0.5*(_.det_cov_z+(aa+uu)/_.s2);
     355      _.logLH
     356    };
     357    ////////////////////////////////////////////////////////////////////////////
     358    VMatrix Draw.U_cond_Z(Real unused)
     359    ////////////////////////////////////////////////////////////////////////////
     360    {
     361      VMatrix u01 = Gaussian(_.q_p,1,0,_.sigma);
     362      _.u + CholeskiSolve(_.L_G, u01,"LtP")
     363    };
     364    ////////////////////////////////////////////////////////////////////////////
     365    Real LogLH.Z_cond_U(Polyn ar, Polyn ma, VMatrix z0, VMatrix a0)
     366    ////////////////////////////////////////////////////////////////////////////
     367    {
     368      -_highValue
     369    };
     370    ////////////////////////////////////////////////////////////////////////////
     371    VMatrix Get.Z0(VMatrix u)
     372    ////////////////////////////////////////////////////////////////////////////
     373    {
     374      Sub(u, 1, 1, _.p, 1)
     375    };
     376    ////////////////////////////////////////////////////////////////////////////
     377    VMatrix Get.A0(VMatrix u)
     378    ////////////////////////////////////////////////////////////////////////////
     379    {
     380      Sub(u, _.p+1, 1, _.q, 1)
     381    };
     382    Real _highValue = 10000*sigma;
     383    Polyn   _.ar;
     384    Polyn   _.ma;
     385    VMatrix _.z     = z_;
     386    VMatrix _.z0    = Rand(0,0,_highValue,_highValue);
     387    VMatrix _.a0    = Rand(0,0,_highValue,_highValue);
     388    Real    _.sigma = sigma;
     389    Real    _.s2    = _.sigma*_.sigma;
     390    Real    _.m     = m;
     391    Real    _.p     = p;
     392    Real    _.q     = q;
     393    Real    _.q_p   = _.q+_.p;
     394    Real    _.mxpq  = Max(_.p,_.q);
     395    VMatrix _.a     = z_*_highValue;
     396    Real _.det_cov_z = _highValue;
     397    Real _.det_cov_u = _highValue;
     398    Real _.det_cov_G = _highValue;
     399    Real _.logLH = _highValue
     400    ]];
     401
     402    If(!p & !q, neutral, {
     403
     404    Polyn   _.ar;
     405    Polyn   _.ma;
    315406    VMatrix _.z     = z_;
    316407    Real    _.sigma = sigma;
    … …  
    322413    Real    _.mxpq  = Max(_.p,_.q);
    323414
    324     ////////////////////////////////////////////////////////////////////////////
    325     VMatrix solve.LLt(VMatrix L, VMatrix b)
    326     ////////////////////////////////////////////////////////////////////////////
    327     {
    328       CholeskiSolve(L,b,"PtLLtP")
    329     };
    330     ////////////////////////////////////////////////////////////////////////////
    331     VMatrix solve.L(VMatrix L, VMatrix b)
    332     ////////////////////////////////////////////////////////////////////////////
    333     {
    334       CholeskiSolve(L,b,"PtL")
    335     };
    336     ////////////////////////////////////////////////////////////////////////////
    337     VMatrix solve.Lt(VMatrix L, VMatrix b)
    338     ////////////////////////////////////////////////////////////////////////////
    339     {
    340       CholeskiSolve(L,b,"LtP")
    341     };
    342415
    343416  //Real _trc.1 = { WriteLn("TRACE Eval.Almagro Creating _.cov_u "+Time); 0 };
    … …  
    348421    //WriteLn("TRACE _.cov_zz = "<<_.cov_zz);
    349422      VMatrix cov_aa = Eye(_.q);
    350       Polyn psi = ChangeBF(Expand(_.ma/_.ar,Max(_.p,_.q)))*B^Max(0,_.q-_.p);
    351       VMatrix cov_za = If(_.q>=_.p, Pol2VMat(psi, _.q, _.p),
    352                             Zeros(_.q, _.p-_.q) | Pol2VMat(psi, _.q, _.q));
    353       VMatrix cov_az = Tra(cov_za);
     423
     424      Polyn psi = Expand(_.ma/_.ar,Max(_.p,_.q));
     425
     426      Polyn psi_F = ChangeBF(psi)*B^Max(0,_.q-_.p);
     427      VMatrix cov_az = If(_.q>=_.p, Pol2VMat(psi_F, _.q, _.p),
     428                            Zeros(_.q, _.p-_.q) | Pol2VMat(psi_F, _.q, _.q));
     429      VMatrix cov_za = Tra(cov_az);
     430/*
     431      VMatrix cov_za_ = Mat2VMat(ForMat(_.p,_.q,Real(Real a, Real b)
     432      {
     433        Real k = _.p + 1 - a;
     434        Real i = _.q + 1 - b;
     435        Coef(psi,i-k)
     436      }));
     437     VMatrix cov_za_dif = cov_za-cov_za_;
     438     Real cov_za.ok = 1-VMatMax(Abs(cov_za_dif));
     439*/
    354440    //WriteLn("TRACE _.cov_za = "<<_.cov_za);
    355441
    356       (cov_zz | cov_az ) <<
    357       (cov_za | cov_aa )
     442      (cov_zz | cov_za ) <<
     443      (cov_az | cov_aa )
    358444    },
    359445    _.p & !_.q,
    … …  
    369455  //WriteLn("_.cov_u = "<<_.cov_u);
    370456    VMatrix _.L_u  = ARMAProcess::CholeskiFactorCov(_.cov_u);
     457    If(!VRows(_.L_u), {
     458    WriteLn("[ARMAProcess::Eval.Almagro] Cannot create _.L_u with \n"+
     459      "AR="<<_.ar+"\n"+
     460      "MA="<<_.ma+"\n");
     461    nonStationary
     462    },{
    371463  //Real _trc.3 = { WriteLn("TRACE Eval.Almagro Creating _.covi_u "+Time); 0 };
    372464    VMatrix _.covi_u =
    373465    {
    374       aux = solve.LLt(_.L_u,Eye(_.q_p));
     466      aux = CholeskiSolve(_.L_u,Eye(_.q_p),"PtLLtP");
    375467      Drop(ARMAProcess::AvoidSymRoundErr(aux),drop)
    376468    };
    … …  
    396488  //Real _trc.6 = { WriteLn("TRACE Eval.Almagro Creating _.L_G "+Time); 0 };
    397489    VMatrix _.L_G   = ARMAProcess::CholeskiFactorCov(_.cov_G);
     490
     491    If(!VRows(_.L_G), {
     492    WriteLn("[ARMAProcess::Eval.Almagro] Cannot create _.L_G with \n"+
     493      "AR="<<_.ar+"\n"+
     494      "MA="<<_.ma+"\n");
     495    nonStationary
     496    },{ [[
     497
     498    ////////////////////////////////////////////////////////////////////////////
     499    Real LogLH.Z(VMatrix z)
     500    ////////////////////////////////////////////////////////////////////////////
     501    {
     502      VMatrix _.z := z;
     503      VMatrix _.pi_z := DifEq(_.ar/_.ma,_.z);
     504      VMatrix _.G_pi_z := Tra(Tra(_.pi_z)*_.G);
     505    //Real _trc.8 = { WriteLn("TRACE Eval.Almagro Creating various "+Time); 0 };
     506      VMatrix _.u := CholeskiSolve(_.L_G, _.G_pi_z,"PtLLtP");
     507      VMatrix _.z0 := Get.Z0(_.u);
     508      VMatrix _.a0 := Get.A0(_.u);
     509      VMatrix _.v := CholeskiSolve(_.L_u, _.u,"PtL");
     510      VMatrix _.a := DifEq(_.ar/_.ma, _.z, _.z0, _.a0);
     511
     512      Real aa = MatDat(VMat2Mat( Tra(_.a)*_.a), 1,1);     
     513      Real uu = MatDat(VMat2Mat( Tra(_.u)*_.covi_u*_.u),1,1);
     514      Real _.logLH := _0.5_m_log2_PI_s2 - 0.5*(_.det_cov_z+(aa+uu)/_.s2);
     515      _.logLH
     516    };
     517    ////////////////////////////////////////////////////////////////////////////
     518    VMatrix Draw.U_cond_Z(Real unused)
     519    ////////////////////////////////////////////////////////////////////////////
     520    {
     521      VMatrix u01 = Gaussian(_.q_p,1,0,_.sigma);
     522      _.u + CholeskiSolve(_.L_G, u01,"LtP")
     523    };
     524
     525    ////////////////////////////////////////////////////////////////////////////
     526    Real LogLH.Z_cond_U(Polyn ar, Polyn ma, VMatrix z0, VMatrix a0)
     527    ////////////////////////////////////////////////////////////////////////////
     528    {
     529      VMatrix a = DifEq(ar/ma,_.z, z0, a0);
     530      _0.5_m_log2_PI_s2 - 0.5*_.m*VMatMoment(a,2)/_.s2
     531    };
     532    ////////////////////////////////////////////////////////////////////////////
     533    VMatrix Get.Z0(VMatrix u)
     534    ////////////////////////////////////////////////////////////////////////////
     535    {
     536      Sub(u, 1, 1, _.p, 1)
     537    };
     538    ////////////////////////////////////////////////////////////////////////////
     539    VMatrix Get.A0(VMatrix u)
     540    ////////////////////////////////////////////////////////////////////////////
     541    {
     542      Sub(u, _.p+1, 1, _.q, 1)
     543    };
     544
     545    Polyn   _.ar;
     546    Polyn   _.ma;
     547    VMatrix _.z;
     548    Real    _.sigma;
     549    Real    _.s2;
     550    Real    _.m;
     551    Real    _.p;
     552    Real    _.q;
     553    Real    _.q_p;
     554    Real    _.mxpq;
     555    VMatrix _.cov_u;
     556    VMatrix _.L_u;
     557    VMatrix _.covi_u;
     558    VMatrix _.L_u_sp;
     559    Real _.det_cov_u;
     560    VMatrix _.H1;
     561    VMatrix _.H2;
     562    VMatrix _.H;
     563    VMatrix _.G;
     564    VMatrix _.GtG;
     565    VMatrix _.cov_G;
     566    VMatrix _.L_G;
     567
     568
    398569  //Real _trc.7 = { WriteLn("TRACE Eval.Almagro Creating _.L_G_sp "+Time); 0 };
    399570    VMatrix _.L_G_sp = Drop(Convert(_.L_G,"Cholmod.R.Sparse"),drop);
    … …  
    401572    Real _.det_cov_z = _.det_cov_u + _.det_cov_G;
    402573    Real _0.5_m_log2_PI_s2 = -0.5*_.m*Log(2*PI*_.s2);
    403 
    404     ////////////////////////////////////////////////////////////////////////////
    405     VMatrix Get.Z0(VMatrix u)
    406     ////////////////////////////////////////////////////////////////////////////
    407     {
    408       Sub(u, 1, 1, _.p, 1)
    409     };
    410     ////////////////////////////////////////////////////////////////////////////
    411     VMatrix Get.A0(VMatrix u)
    412     ////////////////////////////////////////////////////////////////////////////
    413     {
    414       Sub(u, _.p+1, 1, _.q, 1)
    415     };
    416 
    417 
    418574    VMatrix _.pi_z = Rand(0,0,0,0);
    419575    VMatrix _.G_pi_z = Rand(0,0,0,0);
    … …  
    424580    VMatrix _.a = Rand(0,0,0,0);
    425581
    426     Real _.logLH = ?;
    427 
    428     ////////////////////////////////////////////////////////////////////////////
    429     Real LogLH.Z(VMatrix z)
    430     ////////////////////////////////////////////////////////////////////////////
    431     {
    432       VMatrix _.z := z;
    433       VMatrix _.pi_z := DifEq(_.ar/_.ma,_.z);
    434       VMatrix _.G_pi_z := Tra(Tra(_.pi_z)*_.G);
    435     //Real _trc.8 = { WriteLn("TRACE Eval.Almagro Creating various "+Time); 0 };
    436       VMatrix _.u := solve.LLt(_.L_G, _.G_pi_z);
    437       VMatrix _.z0 := Get.Z0(_.u);
    438       VMatrix _.a0 := Get.A0(_.u);
    439       VMatrix _.v := solve.L(_.L_u, _.u);
    440       VMatrix _.a := DifEq(_.ar/_.ma, _.z, _.z0, _.a0);
    441 
    442       Real aa = MatDat(VMat2Mat( Tra(_.a)*_.a), 1,1);     
    443       Real uu = MatDat(VMat2Mat( Tra(_.u)*_.covi_u*_.u),1,1);
    444       Real _.logLH := _0.5_m_log2_PI_s2 - 0.5*(_.det_cov_z+(aa+uu)/_.s2);
    445       _.logLH
    446     };
    447   //Real _trc.9 = { WriteLn("TRACE Eval.Almagro End "+Time); 0 };
    448 
    449     Real _logLH = LogLH.Z(_.z);
    450 
    451     ////////////////////////////////////////////////////////////////////////////
    452     VMatrix Draw.U_cond_Z(Real unused)
    453     ////////////////////////////////////////////////////////////////////////////
    454     {
    455       VMatrix u01 = Gaussian(_.q_p,1,0,_.sigma);
    456       _.u + solve.Lt(_.L_G, u01)
    457     };
    458 
    459582    VMatrix _R0100 = Rand(0,  1,0,0);
    460583    VMatrix _Rq100 = Rand(_.q,1,0,0);
    461584
    462     ////////////////////////////////////////////////////////////////////////////
    463     Real LogLH.Z_cond_U(Polyn ar, Polyn ma, VMatrix z0, VMatrix a0)
    464     ////////////////////////////////////////////////////////////////////////////
    465     {
    466       VMatrix a = DifEq(ar/ma,_.z, z0, a0);
    467       _0.5_m_log2_PI_s2 - 0.5*_.m*VMatMoment(a,2)/_.s2
    468     }
    469   ]])}
     585    Real _.logLH = ?;
     586  //Real _trc.9 = { WriteLn("TRACE Eval.Almagro End "+Time); 0 };
     587    Real _logLH = LogLH.Z(_.z)
     588
     589
     590  ]]}) }) }) }
     591
    470592]];
    471593
Note: See TracChangeset for help on using the changeset viewer.