Changeset 901aaa9 in tol

Timestamp:

Apr 13, 2010, 5:32:52 PM (17 years ago)

Author:

Víctor de Buen <vdebuen@…>

Branches:

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

Children:

e613e67

Parents:

a30612b

Message:

Fixed bug in single try

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

File:

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

tol_tests/tol/ARIMA/mcmc/mtrymh.tol

@@ -159 +159 @@
     VMatrix y = Q_draw(_.x);
     Real pi_logDens_y = pi_logDens(y);
+    Real _.pi_logDens_x := pi_logDens(_.x);
     Real _.r := Min(1, If(_.simetricQ,
     {
@@ -167 +168 @@
           (_.pi_logDens_x + Q_logDens(_.x,   y)) )
     }));
-    Real If(accept_reject(y), _.pi_logDens_x := pi_logDens_y);
+    Real accept_reject(y);
     _.x
   };
@@ -225 +226 @@
       Y[j]
     };
-    VMatrix y = draw_y(0);
-
+    VMatrix y = If(_.k==1,Y[1],draw_y(0));
+  //WriteLn("TRACE draw pi_matchesDomain(y)="<<pi_matchesDomain(y));
     Set X = For(1,_.k-1, VMatrix X_Q_draw(Real j) { Q_draw(y) }) << [[_.x]];
-
     Real calc_denominator_Xjy(Real void)
     {
-    //WriteLn("TRACE draw pi_matchesDomain(y)="<<pi_matchesDomain(y));
       VMatrix lwXjy := Mat2VMat(SetCol(For(1,_.k, Real(Real j)
       {
@@ -242 +241 @@
       VMatrix gXjy = Exp(lwXjy-lwXjy_max);
       Real gXjy_sum := VMatSum(gXjy);
+    //WriteLn("TRACE draw gXjy_sum=\n"<<gXjy_sum);
       True
     };
     Real calc_denominator_Xjy(0);
 
-    //WriteLn("TRACE draw gXjy_sum=\n"<<gXjy_sum);
     Real _.r := Min(1,Exp(lwYjx_max-lwXjy_max)*(gYjx_sum/gXjy_sum));
   //WriteLn("TRACE draw _.r=\n"<<_.r);
     Real accept_reject(y);
-    y
-  })}
+    _.x
+  }) }
 };
 
@@ -616 +615 @@
   ////////////////////////////////////////////////////////////////////////////
   {
+/* * /
+    Real pi_logDens_x = _.pi_logDens_x;
     VMatrix x = _.x;
-    Real pi_logDens_x = _.pi_logDens_x;
-/*
     Matrix Y = Group("ConcatRows", For(1,maxIter,Matrix(Real iter)
     {
@@ -637 +636 @@
       Copy(VMat2Mat(Tra(x)))
     }));
-*/
+/* */
     Matrix pi_logDens_y = Rand(1,maxIter,0,0);
     Matrix Y = Group("ConcatRows", For(1,maxIter,Matrix(Real iter)
@@ -643 +642 @@
       VMatrix y = arima2col(RandCopyARIMA(_.arima));
       Real PutMatDat(pi_logDens_y, 1, iter, pi_logDens(y));
-      Copy(VMat2Mat(Tra(x)))
+      Copy(VMat2Mat(Tra(y)))
     }));
     Real pi_logDens_max = MatMax(pi_logDens_y);
@@ -649 +648 @@
     Matrix g = w * (1/MatSum(w));
     VMatrix _.x := Mat2VMat(Tra(g*Y));
-    Real _.pi_logDens_x := pi_logDens_x;
+    Real _.pi_logDens_x := pi_logDens(_.x);
     Real doAfterDraw(0);
+/* */
     Y
   }
@@ -656 +656 @@
 };
 
+Real PutRandomSeed(0);
+Real rndSeed_iv = GetRandomSeed(0);
+WriteLn("Current Random Seed for initial values= "<<rndSeed_iv);
 
 @MTryMetropolis_ARIMA mtm_arima = @MTryMetropolis_ARIMA::New(
@@ -671 +674 @@
 /* */
 
-Real mtm_arima::_.k := Ceil((mtm_arima::_.n));
+Real mtm_arima::_.k := 9; //Ceil(Sqrt(mtm_arima::_.n));
 
 Matrix mcmc = Rand(0,mtm_arima::_.n,0,0);
@@ -677 +680 @@
 Real mtm_arima::_.s := 0.01;
 
-//Matrix m2avr = mtm_arima::move_to_average(Round(500*Sqrt(mtm_arima::_.n)));
+/*
+Real m2avr.numSim = 500;
+Real m2avr.t0 = Copy(Time);
+Matrix m2avr = mtm_arima::move_to_average(m2avr.numSim);
+Real m2avr.time = (Copy(Time) - m2avr.t0)/m2avr.numSim;
+WriteLn("Time by AVR estimation "<<m2avr.time);
+*/
 
 Real numSim = 0;
@@ -683 +692 @@
 Real append_draw(VMatrix draw)
 {
-  If(!(numSim%10), WriteLn("
+  If(!(numSim%100), WriteLn("
 ##############################################################################
 Simulation "<<numSim+ " _.acceptRatio="<<mtm_arima::_.acceptRatio));
@@ -699 +708 @@
 
 /* */
-Real mtm_arima::_.sFactor := 1.0;
+Real mtm_arima::_.sFactor := 1.00;
 Real numSim := 0;
 Real mtm.t0 = Copy(Time);
 Real TolOprProfiler.Enabled := True;
 While(numSim<2500, append_draw(mtm_arima::draw_multiple_try(0)));
-Real TolOprProfiler.Dump("draw_multiple_try.csv");
+Real TolOprProfiler.Dump("draw_multiple_try.k"<<mtm_arima::_.k+".csv");
 Real mtm.time = (Copy(Time) - mtm.t0)/numSim;
 WriteLn("Time by MTM simulation "<<mtm.time);