Changeset 4414

Show
Ignore:
Timestamp:
03/23/12 17:36:25 (6 years ago)
Author:
vdebuen
Message:

Faster calculation

Files:
1 modified

Legend:

Unmodified
Added
Removed
  • tolp/OfficialTolArchiveNetwork/BysMcmc/bsr/gibbs/_probitFilter.tol

    r4210 r4414  
    2929  Matrix _emptyMatrix = Constant(0,0,?); 
    3030  Matrix _.Y = Constant(0,0,?); 
     31  Matrix _.A = Constant(0,0,?); 
     32  Matrix _.B = Constant(0,0,?); 
    3133  Matrix _residuals_average = Constant(0,0,?); 
    3234  Set _.logDensity = Copy(Empty); 
     
    6365    Matrix _zero := Constant(_.m,1,0); 
    6466    Matrix _inf := _zero+1/0; 
     67    Matrix _.A := IfMat(_.Y,_.Y*(   0),Not(_.Y)*(-1/0)); 
     68    Matrix _.B := IfMat(_.Y,_.Y*(+1/0),Not(_.Y)*(   0)); 
    6569 
    6670    Real _.numEval := 0; 
     
    120124  /////////////////////////////////////////////////////////////////////////// 
    121125  { 
    122   //WriteLn("TRACE [@Filter.Probit::eval.Filter.Probit] 1"); 
     126 
     127  //WriteLn("TRACE"+_MID+"eval 1"); 
     128  //WriteLn("TRACE"+_MID+"eval _.externalInfo[3]=\n"<<Matrix VMat2Mat(_.externalInfo[3],1)); 
     129  //WriteLn("TRACE"+_MID+"eval _.externalInfo[4]=\n"<<Matrix SetRow(_.externalInfo[4])); 
    123130    Matrix Xb = VMat2Mat(SubRow(_.externalInfo[3],_.externalInfo[4])); 
    124   //WriteLn("TRACE [@Filter.Probit::eval.Filter.Probit] 2"); 
    125     Matrix _.normTrunc.lower = IfMat(_.Y, - Xb,  -_inf); 
    126   //WriteLn("TRACE [@Filter.Probit::eval.Filter.Probit] 3"); 
    127     Matrix _.normTrunc.upper = IfMat(_.Y, +_inf, - Xb); 
    128   //WriteLn("TRACE [@Filter.Probit::eval.Filter.Probit] 4"); 
    129     Matrix e = RandTruncatedMultNormal(_zero, _one,  
    130       _.normTrunc.lower, _.normTrunc.upper); 
    131   //WriteLn("TRACE [@Filter.Probit::eval.Filter.Probit] 5"); 
     131  //WriteLn("TRACE"+_MID+"eval Xb=\n"<<Matrix Tra(Xb)); 
     132    Matrix e = RandTruncatedMultNormal(Xb, _one, _.A, _.B); 
    132133    Set selUnk = MatQuery::SelectRowsWithValue(Mat2VMat(Not(IsFinite(e))),1); 
    133134    If(Card(selUnk), 
     
    141142      MatDat(e,k,1)) 
    142143    }); 
    143     Matrix y = Xb+e; 
    144 /* 
    145     Matrix pi0 = F01(Xb); 
    146     Matrix pi1 = -Fy0+1; 
    147     Real logDens = MatSum(IfMat(_.Y, pi1, pi0)); 
    148     Set Append(_.logDensity, [[logDens]]); 
    149   //WriteLn("TRACE [@Filter.Probit::eval.Filter.Probit] 6"); 
    150 /* 
    151     Real burnin = (_.externalInfo[1])::_.config::mcmc.burnin; 
    152     Real cacheLength = (_.externalInfo[1])::_.config::mcmc.cacheLength; 
    153     Real _.numEval := _.numEval + 1; 
    154     Real k = _.numEval%cacheLength; 
    155     Real k_ = If(k==0,cacheLength,k); 
    156   //WriteLn("TRACE [@Filter.Probit::eval.Filter.Probit] 7"); 
    157     Matrix _residuals_average:=If(k==1,e,((_residuals_average*(k_-1))+e)*1/k_); 
    158   //WriteLn("TRACE [@Filter.Probit::eval.Filter.Probit] 8"); 
    159     Set If(k==0, Append(_.residuals_cache,IncludeText( 
    160       "Matrix residuals_average_"<<_.numEval+"=Copy(_residuals_average)"))); 
    161   //WriteLn("TRACE [@Filter.Probit::eval.Filter.Probit] 9"); 
    162 */ 
    163     _.Y-y 
     144  //WriteLn("TRACE"+_MID+"eval e=\n"<<Matrix Tra(e)); 
     145    _.Y-e 
    164146  }; 
    165147