Changeset ffb0c31 in tol

Timestamp:

Nov 13, 2009, 7:12:57 PM (17 years ago)

Author:

Víctor de Buen <vdebuen@…>

Branches:

1941-docker-image-to-test, devel, master

Children:

de893de

Parents:

c255c99

Message:

Storing ARMA initial values

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

Location:

tol/stdlib/math/stat/models/bayesian/bysMcmc

Files:

• _config.tol (modified) (2 diffs)
• bsr/gibbs/_arimaBlock.tol (modified) (7 diffs)

tol/stdlib/math/stat/models/bayesian/bysMcmc/_config.tol

@@ -138 +138 @@
   //Arima filter method
   Text bsr.arimaFilter = BysMcmc::Options::Arima.Filter::FastChol;
+  //Enables/disables storement of the Markov chain of initial values 
+  Real bsr.arimaStoreInitialValues = False;
 //===========================================================================
 //CODA and other diagnostic report configuration
@@ -188 +190 @@
   //If true an internal TOL profiling will be executed
   Real do.mcmc.profiler = False;
+ 
+  @BysInf.Report.Config GetReportConfig(Real void)
+  {
+    @BysInf.Report.Config
+    (
+      report.raftery.diag.q,
+      report.raftery.diag.r,
+      report.raftery.diag.s,
+      report.raftery.diag.eps,
+      report.acf.lag,
+      report.histogram.parts,
+      report.kerDens.points
+    )
+  };
 
   ////////////////////////////////////////////////////////////////////////////

tol/stdlib/math/stat/models/bayesian/bysMcmc/bsr/gibbs/_arimaBlock.tol

@@ -34 +34 @@
   VMatrix   _.difFilter      = Rand(0,0,0,0);
   Set       _.segment.free.index = Copy(Empty);
+  VMatrix   _.U0 = Rand(0,1,0,0); 
+  Matrix    _.U0.cache = Rand(0,0,0,0); 
+  Matrix    _.U0.mcmc = Rand(0,0,0,0); 
+  Set       _.U0.names = Copy(Empty);
 
   //--------------------------------------------------------------------------
@@ -86 +90 @@
     Real numDifData = m;
     VMatrix D = Eye(m);
+    VMatrix u0 = Rand(0,0,0,0);
     Set segment->ArimaAuxInfo := 
     [[
-      numPar, m, difDegree, p, q, numDifData, dif, ar, ma, D
+      numPar, m, difDegree, p, q, numDifData, dif, ar, ma, D, u0
     ]];
     Real SetIndexByName(segment->ArimaAuxInfo);
@@ -191 +196 @@
     Set Li_ar_factors = Copy(Empty);
     Set L_ma_factors = Copy(Empty);
-  //WriteLn("TRACE"+_MID+" initializeArimaBlock["<<resBlockNum_+"] 2.2.7 "+Time);
+    VMatrix u0 = Rand(0,0,0,0);
+    WriteLn("TRACE"+_MID+" initializeArimaBlock 3 ");
+    Set If(_.config::bsr.arimaStoreInitialValues, 
+    {
+      WriteLn("TRACE"+_MID+" initializeArimaBlock 4.1 p="<<p);
+      Set z0.names = For(1,p,Text(Real k)
+      {
+        Text name = segment->Name+".InitialValue.Noise."+IntText(k);
+        Eval(name+"=name")
+      });
+      WriteLn("TRACE"+_MID+" initializeArimaBlock 4.2 q="<<q);
+      Set a0.names = For(1,q,Text(Real k)
+      {
+        Text name = segment->Name+".InitialValue.Residual."+IntText(k);
+        Eval(name+"=name")
+      });
+      Set u0.names = z0.names<<a0.names;
+      WriteLn("TRACE"+_MID+" initializeArimaBlock after 4.3 u0.names=\n"<<u0.names);
+      Append(_.U0.names, u0.names, True)
+    });
+    WriteLn("TRACE"+_MID+" initializeArimaBlock 5 ");
     Set segment->ArimaAuxInfo := 
     [[
        numPar, m, difDegree, p, q, numDifData, 
        dif, ar, ma, Li_ar, L_ma ,Li_ar_factors, L_ma_factors,
-       D, equIdxVec
+       D, equIdxVec, u0
     ]];
     Real SetIndexByName(segment->ArimaAuxInfo);
@@ -538 +563 @@
   //WriteLn("TRACE"+_MID+" draw.scalar [iter "<<numSim+"] 1 ");
     Matrix noise_ = VMat2Mat(noise);
+    Real oldResBlockNum = -1;
+    VMatrix If(_.config::bsr.arimaStoreInitialValues, 
+    {
+      VMatrix _.U0 := Rand(0,1,0,0)
+    }); 
     Set arma.draw.all = 
     EvalSet(_.segment.free.index, Real arma.arms.draw(Real paramNum)
@@ -584 +614 @@
       NameBlock almagro = ARMAProcess::Eval.Almagro(ar,ma,z_,sigma);
 
-      VMatrix u0 = almagro::Draw.U_cond_Z(0);
+      VMatrix u0 = If(oldResBlockNum == pd->resBlockNum,
+      {
+        aux::u0
+      },
+      {
+        aux::u0 := almagro::Draw.U_cond_Z(0);
+        If(_.config::bsr.arimaStoreInitialValues, 
+        {
+          _.U0 := _.U0 << aux::u0
+        });
+        aux::u0  
+      });  
       VMatrix z0 = almagro::Get.Z0(u0);
       VMatrix a0 = almagro::Get.A0(u0);
@@ -609 +650 @@
     //WriteLn("TRACE"+_MID+" draw.scalar pd->coef="<<pd->coef);
       pd->coef
+    });
+
+    Matrix If(_.config::bsr.arimaStoreInitialValues,
+    {
+      Case(
+      numSim==1,
+      {
+        _.U0.mcmc := Rand(0,VRows(_.U0),0,0);
+        _.U0.cache := VMat2Mat(_.U0,True)
+      },
+      numSim%100,
+      {
+        _.U0.cache := _.U0.cache << VMat2Mat(_.U0,True)
+      },
+      1==1,
+      {
+        _.U0.mcmc := _.U0.mcmc << _.U0.cache;
+        _.U0.cache := Rand(0,VRows(_.U0),0,0)
+      })
     });
   //Real show.arma(0);
@@ -803 +863 @@
     Set _.arimaStructBlock = 
     {
+      BysMcmc::@Config _.config := Copy(config);
+
     //WriteLn("TRACESelecting noise with known wariance");
       Select(NoiseDistrib, Real(@BSR.NoiseDistrib segment)