Changeset f9c7f71 in tol

Timestamp:

Apr 14, 2010, 3:51:56 PM (17 years ago)

Author:

Víctor de Buen <vdebuen@…>

Branches:

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

Children:

823054d

Parents:

fdb5742

Message:

La llamada a Estimate es opcional
Se reporta la comparación de verosimilitudes del ARIMA generado real y el del valor medio de la MCMC

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

File:

• tol_tests/tol/ARIMA/mcmc/mtrymh.mcmc.tol (modified) (9 diffs)

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

@@ -10 +10 @@
 Matrix mcmc = Rand(0,mtm_arima::_.n,0,0);
 
-Matrix L = {
-  Matrix COV = modEst.0::Correlations::COV;
-  Real n = mtm_arima::_.n;
-  Real r = Rows(COV);
-  Real k = r-n+1;
-  Choleski(Sub(COV,k,k,n,n))
-};
-Matrix Li = LTInverse(L);
-VMatrix mtm_arima::_.L := Mat2VMat(L);
-VMatrix mtm_arima::_.Li := Mat2VMat(Li);
-Real mtm_arima::_.s := 1;
-
+Real If(!ObjectExist("Matrix","modEst.0::Correlations::COV"),
+{
+  Real mtm_arima::_.s := 0.001
+},
+{
+  Matrix L = {
+    Matrix COV = modEst.0::Correlations::COV;
+    Real n = mtm_arima::_.n;
+    Real r = Rows(COV);
+    Real k = r-n+1;
+    Choleski(Sub(COV,k,k,n,n))
+  };
+  Matrix Li = LTInverse(L);
+  VMatrix mtm_arima::_.L := Mat2VMat(L);
+  VMatrix mtm_arima::_.Li := Mat2VMat(Li);
+  Real mtm_arima::_.s := 1
+});
+  
 
 Real numSim = 0;
@@ -43 +49 @@
 
 /* */
-Real mtm_arima::_.sFactor := 1.00;
+Real mtm_arima::_.sFactor := 1.01;
 Real numSim := 0;
 Real mtm.t0 = Copy(Time);
@@ -62 +68 @@
 Real raftery.rate = ?;
 Real raftery.speed = ?;
+Real burnin = 500;
 
 Matrix cmp = mcmc | Group("ConcatColumns", For(1,Rows(PARAM),Matrix(Real j) 
@@ -78 +85 @@
   ]]});
   Real raftery.burnin := MatMax(SubCol(raftery,[[1]]));
+  Real burnin := Min(raftery.burnin,500);
   Real raftery.thinning := MatMax(SubCol(raftery,[[4]]));
   Real raftery.rate := (numSim-raftery.burnin)/(numSim*raftery.thinning);
@@ -83 +91 @@
 ]]});
   
+Matrix mcmc.sub = SubRow(mcmc, Range(burnin+1,Rows(mcmc),1));
+Matrix x.avr = Tra(Constant(1,Rows(mcmc.sub),1/Rows(mcmc.sub))*mcmc.sub);
+Set arima.avr = mtm_arima::col2arima(Mat2VMat(x.avr));
+Set ev.avr = ARIMALevinsonEval(arima.avr, Tra(SerMat(z)));
+Serie e;
+Serie residuals.avr  = 
+{
+  Matrix res = ev.avr[4];
+  Date   ini = Succ(Last(z),Diario,-Rows(ev.avr[4])+1);
+  If(!Rows(res), 0*z, MatSerSet(Tra(res),Diario,ini)[1])
+};
+Real logLH.real = (SimulArima::EvalRealARIMA::ev)[12];
+Real logLH.avr = (ev.avr)[12];
+Real prob.refuse = Min(1,Max(0,1-Exp(logLH.avr-logLH.real)));
+
+Real Q_draw_by_Pi_draw = mtm_arima::_.num_Q_draw/mtm_arima::_.num_draw;
+Real num_logWeight_by_Pi_draw = mtm_arima::_.num_logWeight/mtm_arima::_.num_draw;
+
 Text report = ""
   <<rndSeed+"\t"
@@ -89 +115 @@
   <<mtm_arima::_.n+"\t"
   <<mtm_arima::_.k+"\t"
+  <<Q_draw_by_Pi_draw+"\t"
+  <<num_logWeight_by_Pi_draw+"\t"
   <<arima.time+"\t"
   <<mtm.time+"\t"
@@ -97 +125 @@
   <<raftery.rate+"\t"
   <<raftery.speed+"\t"
+  <<logLH.real+"\t"
+  <<logLH.avr+"\t"
+  <<prob.refuse+"\t"
   <<numErr+"\t"
   <<numWar+"\n";
@@ -105 +136 @@
   +"mtm_arima::_.n\t"
   +"mtm_arima::_.k\t"
+  +"Q_draw_by_Pi_draw\t"
+  +"num_logWeight_by_Pi_draw\t"
   +"arima.time\t"
   +"mtm.time\t"
@@ -112 +145 @@
   +"raftery.thinning\t"
   +"raftery.rate\t"
-  +"raftery.speed\t"
+  +"logLH.real\t"
+  +"logLH.avr\t"
+  +"prob.refuse\t"
   +"numErr\t"
   +"numWar\n"));