| 1 | /* -*- mode: C++ -*- */
|
|---|
| 2 | //////////////////////////////////////////////////////////////////////////////
|
|---|
| 3 | // FILE : test.tol
|
|---|
| 4 | //////////////////////////////////////////////////////////////////////////////
|
|---|
| 5 |
|
|---|
| 6 |
|
|---|
| 7 | Text email_ = "vdebuen@tol-project.org";
|
|---|
| 8 | Text link_ = "";
|
|---|
| 9 | Text summary_ =
|
|---|
| 10 | "test_0001 multivariate linear model with known covariance matrix" ;
|
|---|
| 11 |
|
|---|
| 12 | Real numErr0 = Copy(NError);
|
|---|
| 13 |
|
|---|
| 14 | #Embed "../gen_multivariant_linear_regression.tol";
|
|---|
| 15 |
|
|---|
| 16 |
|
|---|
| 17 | BysVecLinReg::@Multivariate MV = [[
|
|---|
| 18 | Text _.desc = "Modelo lineal multivariante para ventas de Fanta en Turquía",
|
|---|
| 19 | Set _.outputNames = For(1,d,Text(Real i) { "Output_"<<i }),
|
|---|
| 20 | Set _.inputNames = For(1,m,Text(Real i) { "Input_"<<i })
|
|---|
| 21 | ]];
|
|---|
| 22 |
|
|---|
| 23 | Real MV::set.outputInput(Y,X);
|
|---|
| 24 | Real MV::change.V(V);
|
|---|
| 25 |
|
|---|
| 26 | VMatrix S = KroneckerProduct(V,MtMSqr(MV::_.X.Li));
|
|---|
| 27 | VMatrix L = Tra(MV::_.IdxXLi) * MV::_.VxIm.Le;
|
|---|
| 28 | VMatrix LLt = MMtSqr(L);
|
|---|
| 29 |
|
|---|
| 30 | Real ok_.X.Li = VMatMax(Abs((MV::_.X.Li * MV::_.X.Le+m*N)-m*N - Eye(m)))==0;
|
|---|
| 31 | Real ok_.IdxXL = VMatMax(Abs((Tra(MV::_.IdxXL)*S*MV::_.IdxXL - MV::_.VxIm+m*d)-m*d))==0;
|
|---|
| 32 | Real ok_.IdxXLi = VMatMax(Abs((MV::_.IdxXLi * MV::_.IdxXL+m*N)-m*N - Eye(m*d)))==0;
|
|---|
| 33 | Real okLLt = VMatMax(Abs((LLt-S+m*N)-m*N ))==0;
|
|---|
| 34 |
|
|---|
| 35 | VMatrix mcmc = Group("ConcatRows", For(1,sampleLength,VMatrix(Real numSim)
|
|---|
| 36 | {
|
|---|
| 37 | Tra(MV::gibbs.draw.vec.A(0))
|
|---|
| 38 | }));
|
|---|
| 39 |
|
|---|
| 40 | VMatrix vec.A_ = Tra(Constant(1,sampleLength,1/sampleLength)*mcmc);
|
|---|
| 41 | VMatrix vec.A = BysVecLinReg::vmatrix2vec(A);
|
|---|
| 42 | VMatrix cmp.A = vec.A | vec.A_;
|
|---|
| 43 |
|
|---|
| 44 | VMatrix A_ = BysVecLinReg::vec2vmatrix(vec.A_,A);
|
|---|
| 45 | VMatrix E_ = Y - A_*X;
|
|---|
| 46 | VMatrix E01_ = CholeskiSolve(V.L,E_,"PtLLtP");
|
|---|
| 47 | VMatrix cmpE01 = BysVecLinReg::vmatrix2vec(E01) |
|
|---|
| 48 | BysVecLinReg::vmatrix2vec(E01_);
|
|---|
| 49 |
|
|---|
| 50 | Real stdErr = VMatAvr(E01 ^2);
|
|---|
| 51 | Real stdErr_ = VMatAvr(E01_^2);
|
|---|
| 52 |
|
|---|
| 53 | Real numErr1 = Copy(NError);
|
|---|
| 54 |
|
|---|
| 55 | //Extented information about test behaviour
|
|---|
| 56 | Set partialResults_ =
|
|---|
| 57 | [[
|
|---|
| 58 | numErr0, numErr1, stdErr, stdErr_,
|
|---|
| 59 | ok_.X.Li, ok_.IdxXL, ok_.IdxXLi, okLLt
|
|---|
| 60 |
|
|---|
| 61 | ]];
|
|---|
| 62 |
|
|---|
| 63 | Real quality_ = And(
|
|---|
| 64 | numErr1==numErr0,
|
|---|
| 65 | ok_.X.Li, ok_.IdxXL, ok_.IdxXLi, okLLt)*
|
|---|
| 66 | Min(1,Exp((stdErr-stdErr_)*N*d));
|
|---|
| 67 |
|
|---|
| 68 | //Return the results
|
|---|
| 69 | Set resultStr_ = @strTestStatus(summary_, link_, quality_,
|
|---|
| 70 | "Partial results = "<<partialResults_,
|
|---|
| 71 | "NO DBServerType", "NO DBServerHost", "NO DBAlias",
|
|---|
| 72 | email_);
|
|---|
| 73 | WriteLn(""<<resultStr_);
|
|---|
| 74 | resultStr_
|
|---|
| 75 |
|
|---|
| 76 | /* */
|
|---|
| 77 |
|
|---|
| 78 |
|
|---|
| 79 | /* */
|
|---|
