Changeset 8ffe39c in tol

Timestamp:

Jan 5, 2012, 11:24:50 AM (15 years ago)

Author:

Víctor de Buen <vdebuen@…>

Branches:

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

Children:

16db681

Parents:

7b5e14e

Message:

Refs #1387

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

Location:

tol_tests/tol/BSR/test_08

Files:

• fit_beta.tol (added)
• fit_beta_01.gif (added)
• fit_beta_02.gif (added)
• fit_beta_03.gif (added)
• fit_beta_04.gif (added)
• test.tol (modified) (9 diffs)

tol_tests/tol/BSR/test_08/test.tol

@@ -26 +26 @@
 Matrix zero = Rand(m,1,0,0);
 Matrix one = zero+1;
-//Original supplied data 
-Matrix X1 = Rand(m, 5, -.3, .3);
-//Deived inputs
-Matrix X2 = SubCol(X1,[[1]]) $* SubCol(X1,[[2]]);
-//Constant input
-Matrix X3 = Rand(Rows(X1),1,1,1);
 //Full input matrix
-Matrix X  = X1|X2|X3;
-//Real parameters of linear block
+Matrix X  = {
+  Matrix X1 = Round(Rand(m, 5, -0.5, 1.5));
+  //Matrix X2 = SubCol(X1,[[1]]) $* SubCol(X1,[[2]]);
+  Matrix X2 = Round(Rand(m, 1, -0.5, 2.5));
+  //Constant input
+  Matrix X3 = Rand(m,1,1,1);
+  Matrix X_ = X1|X2|X3
+};
+
+  //Real parameters of linear block
 Matrix b  = {Col(
-  Param_1=Rand(-2,0),
-  Param_2=Rand(-2,0),
+  Param_1=Rand(-1,0),
+  Param_2=Rand( 0,1),
   Param_3=Rand(-2,0),
-  Param_4=Rand( 0,1),
+  Param_4=Rand( 0,0.2),
   Param_5=Rand(Param_4,2),
   Param_6=Rand(Param_5,3),
@@ -46 +48 @@
 Matrix e  = Gaussian(Rows(X),1,0,1);
 //Real filter matrix
-Matrix F = X*b;
+Matrix Xb = X*b;
 //Latent linear output
-Matrix Z = F+e;
+Matrix z = Xb+e;
 //Real output probabilities
-Matrix Yp = F01(Z);
+Matrix Yp = F01(Xb);
 //Simulated probit model
-Matrix Y  = GE(Z,zero);
+Matrix Y  = GE(z,zero);
+Real Y_avr = MatAvr(Y);
+Matrix Xb_e = Xb | e;
+Matrix Yp_Ype = Yp | F01(Xb+e);
+
+Real success  = MatAvr(EQ(Y,GE(Yp,    0.5)));
+
+/* */
 
 Set mh.estim = Probit(Y, X);
 Matrix b.mh = mh.estim[1];
-
+Matrix Xb.mh = X*b.mh;
+Matrix Yp.mh = F01(Xb.mh);
 //////////////////////////////////////////////////////////////////////////////
 //These is the only information that will be known in real world problems
 //////////////////////////////////////////////////////////////////////////////
 //Data matrix containing full output and input
-Matrix data.mat = Y|X1;
+Matrix data.mat = Y|X;
 //Set of names of columns in data.mat
-Set data.name = [["Output"]]<<For(1,Columns(X1),Text(Real j)
+Set data.name = [["Output"]]<<For(1,Columns(X),Text(Real j)
 { "Input_"+IntText(j)} );
 
@@ -98 +108 @@
   Real _check_data.orig = SetIndexByName(_.data.orig);
 
-  Set _.data.handler = _.data.orig <<
-  {[[
-    Matrix Input_Prod_1_2 = _.data.orig::Input_1 $* _.data.orig::Input_2;
-    Matrix Constant = Rand(Rows(_.data.orig[1]),1,1,1)
-  ]]};
+  Set _.data.handler = _.data.orig;
 
   Real _check_data.handler = SetIndexByName(_.data.handler);
@@ -119 +125 @@
     //   Variable      initVal   min    max     priAvg   priSig  priSigWeight
     LPI("Input_1",           0, -1/0,     0,         ?,       ?,          ?),
-    LPI("Input_2",           0, -1/0,     0,         ?,       ?,          ?),
+    LPI("Input_2",           0,    0,  +1/0,         ?,       ?,          ?),
     LPI("Input_3",           0, -1/0,     0,         ?,       ?,          ?),
     LPI("Input_4",           0,    0,  +1/0,         ?,       ?,          ?),
     LPI("Input_5",           0,    0,  +1/0,         ?,       ?,          ?),
-    LPI("Input_Prod_1_2",    0,    0,  +1/0,         ?,       ?,          ?),
-    LPI("Constant" ,  _.cteAvg, -1/0,  +1/0,  _.cteAvg,    0.25,       0.90) 
+    LPI("Input_6",           0,    0,  +1/0,         ?,       ?,          ?),
+    LPI("Input_7" ,   _.cteAvg, -1/0,  +1/0,         ?,       ?,          ?) 
   ]]};
   //Input matrix
@@ -140 +146 @@
     //                      lower        <=  upper  
     BysMcmc::@Bsr.OrderRelation.Info("Input_4",       "Input_5"),
-    BysMcmc::@Bsr.OrderRelation.Info("Input_5",       "Input_Prod_1_2")
+    BysMcmc::@Bsr.OrderRelation.Info("Input_5",       "Input_6")
   ]]
 
@@ -151 +157 @@
 //////////////////////////////////////////////////////////////////////////////
 //MCMC dimensions
-  Real mcmc.burnin        =  1000; 
-  Real mcmc.sampleLength  =  4000;
-  Real mcmc.thinning      =     2; //Probit mcmc have high alternance
-  Real mcmc.cacheLength   =   100;
-    Real bsr.iniVal.try.findFeasible = 1;
-    Real bsr.iniVal.try.zero = 2;
-    Real bsr.iniVal.try.givenByUser = 3;
+  Real mcmc.burnin        =   2000; 
+  Real mcmc.sampleLength  =  10000;
+  Real mcmc.thinning      =  1; //Probit mcmc have high alternance
+  Real mcmc.cacheLength   =  If(mcmc.sampleLength>=2000,100,1);
+    Real bsr.iniVal.try.givenByUser = 1;
+    Real bsr.iniVal.try.cnstrMLE = 2;
+    Real bsr.iniVal.try.zero = 3;
+    Real bsr.iniVal.try.findFeasible = 4;
     Real bsr.iniVal.try.cnstrMinNorm = 5;
-    Real bsr.iniVal.try.cnstrMLE = 4;
 //Generic flags
   Real do.resume          =  0;
-  Real do.report          =  True;
+  Real do.report          =  mcmc.sampleLength>=2000;
   Real do.eval            =  False; //Has no sense for probit models
   Real do.linear.effects  =  False  //Has no sense for probit models
@@ -187 +193 @@
                  Rows(b_full)-VColumns(bsr.oneNode.data::_.X),0,0);
 
+Matrix Xb_ = X*b_;
+Matrix Yp_ = F01(Xb_);
+
+Matrix CompareBeta = b | b.mh | b_;
+Matrix CompareFilter = Xb | Xb.mh | Xb_;
+Matrix CompareProb = Yp| Yp.mh | Yp_;
+
+Real success.mh = MatAvr(EQ(Y,GE(Yp.mh, 0.5)));
+Real success_   = MatAvr(EQ(Y,GE(Yp_,   0.5)));
+
+
+Set residuals_cache = ((bsr.oneNode.cycler::_.sampler::_.nonLinBlk)[1])::
+  _.filter::_.residuals_cache;
+
+Set residuals_cache_stats = For(1,Card(residuals_cache),Set(Real k)
+{
+  Matrix ek = residuals_cache[k]; 
+  Real nu = MatAvr(ek);
+  Real sig = Sqrt(MatAvr(ek$*ek));
+  Matrix frq = Frequency(ek,30,-3,3);
+  [[k,nu,sig,frq]]
+});
+Set rc = ExtractByIndex(residuals_cache,Range(10,Card(residuals_cache),1));
+Matrix e_ = SetSum(rc)*1/Card(rc);
+Real nu_ = MatAvr(e_);
+Real sig_ = Sqrt(MatAvr(e_$*e_));
+Matrix frq_ = Frequency(e_,30,-3,3);
+
+
+Real numErr1 = Copy(NError);
+
+Real quality = Min(1,success_/success);
+//Extented information about test behaviour
+Set partialResults_ = 
+[[
+  numErr0, numErr1
+]];
+
+//This is a messure of the success of the test 
+Real quality_ = And(numErr0 == numErr1)*quality;
+
+
+//Return the results 
+Set resultStr_ = @strTestStatus(summary_, link_, quality_,
+                  "Partial results = "<<partialResults_,
+                  "NO DBServerType", "NO DBServerHost", "NO DBAlias",
+                  email_);
+WriteLn(""<<resultStr_);
+resultStr_
+
+/*
 //Output probabilities for given input matrix integrated over all simulated 
 //vectors of parameters (Bayesian inference by exact method)
@@ -197 +254 @@
 Matrix outAprBys = GE(outProbAprBys,one*0.5);
 
-Matrix CompareBeta = b | b_ | b.mh;
-Matrix CompareProb = Yp|outProbExcBys|outProbAprBys;
 Matrix SuccessExcBys = EQ(Y,outExcBys);
 Matrix SuccessAprBys = EQ(Y,outAprBys);
@@ -213 +268 @@
   outProbAprBys);
 
-Real numErr1 = Copy(NError);
-
-Real quality = MatAvr(SuccessExcBys);
-//Extented information about test behaviour
-Set partialResults_ = 
-[[
-  numErr0, numErr1
-]];
-
-//This is a messure of the success of the test 
-Real quality_ = And(numErr0 == numErr1)*quality;
-
-
-//Return the results 
-Set resultStr_ = @strTestStatus(summary_, link_, quality_,
-                  "Partial results = "<<partialResults_,
-                  "NO DBServerType", "NO DBServerHost", "NO DBAlias",
-                  email_);
-WriteLn(""<<resultStr_);
-resultStr_
-
 /* */