nm_misc.cpp
Go to the documentation of this file.
00001 
00002 
00003 
00006 
00007 #define WANT_MATH
00008 
00009 #include "include.h"
00010 
00011 #include "newmatap.h"
00012 
00013 #ifdef use_namespace
00014 namespace NEWMAT {
00015 #endif
00016 
00017 
00018 #ifdef DO_REPORT
00019 #define REPORT { static ExeCounter ExeCount(__LINE__,21); ++ExeCount; }
00020 #else
00021 #define REPORT {}
00022 #endif
00023 
00024 ReturnMatrix Helmert(int n, bool full)
00025 {
00026    REPORT
00027    Tracer et("Helmert ");
00028    if (n <= 0) Throw(ProgramException("Dimension <= 0 "));
00029    Matrix H;
00030    
00031    if (full) H.resize(n,n); else H.resize(n-1,n);
00032    H = 0.0;
00033    for (int i = 1; i < n; ++i)
00034    {
00035       Real f = 1.0 / sqrt((Real)i * (i+1));
00036       H.submatrix(i,i,1,i) = -f; H(i,i+1) = f * i;
00037    }
00038    if (full) { H.row(n) = 1.0 / sqrt((Real)n); }
00039    H.release(); return H.for_return();
00040 } 
00041 
00042 
00043 
00044 // Multiply X by n-1 x n matrix to give n-1 contrasts
00045 // Return a ColumnVector
00046 ReturnMatrix Helmert(const ColumnVector& X, bool full)
00047 {
00048    REPORT
00049    Tracer et("Helmert * CV");
00050    int n = X.nrows();
00051    if (n == 0) Throw(ProgramException("X Vector of length 0", X));
00052    Real sum = 0.0; ColumnVector Y;
00053    if (full) Y.resize(n); else Y.resize(n-1);
00054    for (int i = 1; i < n; ++i)
00055       { sum += X(i); Y(i) = (i * X(i+1) - sum) / sqrt((Real)i * (i+1)); }
00056    if (full) { sum += X(n); Y(n) = sum / sqrt((Real)n); }
00057    Y.release(); return Y.for_return();
00058 } 
00059 
00060 // same as above for X a ColumnVector, length n, element j = 1; otherwise 0
00061 ReturnMatrix Helmert(int n, int j, bool full)
00062 {
00063    REPORT
00064    Tracer et("Helmert:single element ");
00065    if (n <= 0) Throw(ProgramException("X Vector of length <= 0"));
00066    if (j > n || j <= 0)
00067       Throw(ProgramException("Out of range element number "));
00068    ColumnVector Y; if (full) Y.resize(n); else Y.resize(n-1);
00069    Y = 0.0;
00070    if (j > 1) Y(j-1) = sqrt((Real)(j-1) / (Real)j);
00071    for (int i = j; i < n; ++i) Y(i) = - 1.0 / sqrt((Real)i * (i+1));
00072    if (full) Y(n) = 1.0 / sqrt((Real)n);
00073    Y.release(); return Y.for_return();
00074 } 
00075 
00076 ReturnMatrix Helmert_transpose(const ColumnVector& Y, bool full)
00077 {
00078    REPORT
00079    Tracer et("Helmert_transpose * CV ");
00080    int n = Y.nrows(); Real sum;
00081    if (full) sum = Y(n) / sqrt((Real)n); else { sum = 0.0; ++n; }
00082    ColumnVector X(n);
00083    for (int i = n-1; i > 0; --i)
00084    {
00085       Real Yi = Y(i) / sqrt((Real)i * (i+1));
00086       X(i+1) = i * Yi + sum; sum -= Yi;
00087    }
00088    X(1) = sum;
00089    X.release(); return X.for_return();
00090 }
00091 
00092 // same as above but want only j-th element
00093 Real Helmert_transpose(const ColumnVector& Y, int j, bool full)
00094 {
00095    REPORT
00096    Tracer et("Helmert_transpose:single element ");
00097    int n = Y.nrows(); Real sum;
00098    if (full) sum = Y(n) / sqrt((Real)n); else { sum = 0.0; ++n; }
00099    if (j > n || j <= 0) Throw(IndexException(j, Y));
00100    for (int i = n-1; i > 0; --i)
00101    {
00102       Real Yi = Y(i) / sqrt((Real)i * (i+1));
00103       if (i+1 == j) return i * Yi + sum;
00104       sum -= Yi;
00105    }
00106    return sum;
00107 }
00108 
00109 ReturnMatrix Helmert(const Matrix& X, bool full)
00110 {
00111    REPORT
00112    Tracer et("Helmert * Matrix");
00113    int m = X.nrows(); int n = X.ncols();
00114    if (m == 0) Throw(ProgramException("Matrix has 0 rows ", X));
00115    Matrix Y;
00116    if (full) Y.resize(m,n); else Y.resize(m-1, n);
00117    for (int j = 1; j <= n; ++j)
00118    {
00119       ColumnVector CV = X.Column(j);
00120       Y.Column(j) = Helmert(CV, full);
00121    }
00122    Y.release(); return Y.for_return();
00123 }
00124 
00125 ReturnMatrix Helmert_transpose(const Matrix& Y, bool full)
00126 {
00127    REPORT
00128    Tracer et("Helmert_transpose * Matrix ");
00129    int m = Y.nrows(); int n = Y.ncols(); if (!full) ++m;
00130    Matrix X(m, n);
00131    for (int j = 1; j <= n; ++j)
00132    {
00133       ColumnVector CV = Y.Column(j);
00134       X.Column(j) = Helmert_transpose(CV, full);
00135    }
00136    X.release(); return X.for_return();
00137 }
00138 
00139 
00140 
00141 
00142 #ifdef use_namespace
00143 }
00144 #endif
00145 
00146 


kni
Author(s): Martin Günther
autogenerated on Thu Jun 6 2019 21:42:34