本文整理汇总了C++中SymmetricMatrix::i方法的典型用法代码示例。如果您正苦于以下问题:C++ SymmetricMatrix::i方法的具体用法?C++ SymmetricMatrix::i怎么用?C++ SymmetricMatrix::i使用的例子?那么恭喜您, 这里精选的方法代码示例或许可以为您提供帮助。您也可以进一步了解该方法所在类SymmetricMatrix
的用法示例。
在下文中一共展示了SymmetricMatrix::i方法的4个代码示例,这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞,您的评价将有助于系统推荐出更棒的C++代码示例。
示例1: test1
void test1(Real* y, Real* x1, Real* x2, int nobs, int npred)
{
cout << "\n\nTest 1 - traditional, bad\n";
// traditional sum of squares and products method of calculation
// but not adjusting means; maybe subject to round-off error
// make matrix of predictor values with 1s into col 1 of matrix
int npred1 = npred+1; // number of cols including col of ones.
Matrix X(nobs,npred1);
X.Column(1) = 1.0;
// load x1 and x2 into X
// [use << rather than = when loading arrays]
X.Column(2) << x1; X.Column(3) << x2;
// vector of Y values
ColumnVector Y(nobs); Y << y;
// form sum of squares and product matrix
// [use << rather than = for copying Matrix into SymmetricMatrix]
SymmetricMatrix SSQ; SSQ << X.t() * X;
// calculate estimate
// [bracket last two terms to force this multiplication first]
// [ .i() means inverse, but inverse is not explicity calculated]
ColumnVector A = SSQ.i() * (X.t() * Y);
// Get variances of estimates from diagonal elements of inverse of SSQ
// get inverse of SSQ - we need it for finding D
DiagonalMatrix D; D << SSQ.i();
ColumnVector V = D.AsColumn();
// Calculate fitted values and residuals
ColumnVector Fitted = X * A;
ColumnVector Residual = Y - Fitted;
Real ResVar = Residual.SumSquare() / (nobs-npred1);
// Get diagonals of Hat matrix (an expensive way of doing this)
DiagonalMatrix Hat; Hat << X * (X.t() * X).i() * X.t();
// print out answers
cout << "\nEstimates and their standard errors\n\n";
// make vector of standard errors
ColumnVector SE(npred1);
for (int i=1; i<=npred1; i++) SE(i) = sqrt(V(i)*ResVar);
// use concatenation function to form matrix and use matrix print
// to get two columns
cout << setw(11) << setprecision(5) << (A | SE) << endl;
cout << "\nObservations, fitted value, residual value, hat value\n";
// use concatenation again; select only columns 2 to 3 of X
cout << setw(9) << setprecision(3) <<
(X.Columns(2,3) | Y | Fitted | Residual | Hat.AsColumn());
cout << "\n\n";
}
示例2: trymatd
void trymatd()
{
Tracer et("Thirteenth test of Matrix package");
Tracer::PrintTrace();
Matrix X(5,20);
int i,j;
for (j=1;j<=20;j++) X(1,j) = j+1;
for (i=2;i<=5;i++) for (j=1;j<=20; j++) X(i,j) = (long)X(i-1,j) * j % 1001;
SymmetricMatrix S; S << X * X.t();
Matrix SM = X * X.t() - S;
Print(SM);
LowerTriangularMatrix L = Cholesky(S);
Matrix Diff = L*L.t()-S; Clean(Diff, 0.000000001);
Print(Diff);
{
Tracer et1("Stage 1");
LowerTriangularMatrix L1(5);
Matrix Xt = X.t(); Matrix Xt2 = Xt;
QRZT(X,L1);
Diff = L - L1; Clean(Diff,0.000000001); Print(Diff);
UpperTriangularMatrix Ut(5);
QRZ(Xt,Ut);
Diff = L - Ut.t(); Clean(Diff,0.000000001); Print(Diff);
Matrix Y(3,20);
for (j=1;j<=20;j++) Y(1,j) = 22-j;
for (i=2;i<=3;i++) for (j=1;j<=20; j++)
Y(i,j) = (long)Y(i-1,j) * j % 101;
Matrix Yt = Y.t(); Matrix M,Mt; Matrix Y2=Y;
QRZT(X,Y,M); QRZ(Xt,Yt,Mt);
Diff = Xt - X.t(); Clean(Diff,0.000000001); Print(Diff);
Diff = Yt - Y.t(); Clean(Diff,0.000000001); Print(Diff);
Diff = Mt - M.t(); Clean(Diff,0.000000001); Print(Diff);
Diff = Y2 * Xt2 * S.i() - M * L.i();
Clean(Diff,0.000000001); Print(Diff);
}
ColumnVector C1(5);
{
Tracer et1("Stage 2");
X.ReSize(5,5);
for (j=1;j<=5;j++) X(1,j) = j+1;
for (i=2;i<=5;i++) for (j=1;j<=5; j++)
X(i,j) = (long)X(i-1,j) * j % 1001;
for (i=1;i<=5;i++) C1(i) = i*i;
CroutMatrix A = X;
ColumnVector C2 = A.i() * C1; C1 = X.i() * C1;
X = C1 - C2; Clean(X,0.000000001); Print(X);
}
{
Tracer et1("Stage 3");
X.ReSize(7,7);
for (j=1;j<=7;j++) X(1,j) = j+1;
for (i=2;i<=7;i++) for (j=1;j<=7; j++)
X(i,j) = (long)X(i-1,j) * j % 1001;
C1.ReSize(7);
for (i=1;i<=7;i++) C1(i) = i*i;
RowVector R1 = C1.t();
Diff = R1 * X.i() - ( X.t().i() * R1.t() ).t(); Clean(Diff,0.000000001);
Print(Diff);
}
{
Tracer et1("Stage 4");
X.ReSize(5,5);
for (j=1;j<=5;j++) X(1,j) = j+1;
for (i=2;i<=5;i++) for (j=1;j<=5; j++)
X(i,j) = (long)X(i-1,j) * j % 1001;
C1.ReSize(5);
for (i=1;i<=5;i++) C1(i) = i*i;
CroutMatrix A1 = X*X;
ColumnVector C2 = A1.i() * C1; C1 = X.i() * C1; C1 = X.i() * C1;
X = C1 - C2; Clean(X,0.000000001); Print(X);
}
{
Tracer et1("Stage 5");
int n = 40;
SymmetricBandMatrix B(n,2); B = 0.0;
for (i=1; i<=n; i++)
{
B(i,i) = 6;
if (i<=n-1) B(i,i+1) = -4;
if (i<=n-2) B(i,i+2) = 1;
}
B(1,1) = 5; B(n,n) = 5;
SymmetricMatrix A = B;
ColumnVector X(n);
X(1) = 429;
for (i=2;i<=n;i++) X(i) = (long)X(i-1) * 31 % 1001;
X = X / 100000L;
// the matrix B is rather ill-conditioned so the difficulty is getting
// good agreement (we have chosen X very small) may not be surprising;
// maximum element size in B.i() is around 1400
ColumnVector Y1 = A.i() * X;
LowerTriangularMatrix C1 = Cholesky(A);
ColumnVector Y2 = C1.t().i() * (C1.i() * X) - Y1;
Clean(Y2, 0.000000001); Print(Y2);
UpperTriangularMatrix CU = C1.t().i();
//.........这里部分代码省略.........
示例3: test2
void test2(Real* y, Real* x1, Real* x2, int nobs, int npred)
{
cout << "\n\nTest 2 - traditional, OK\n";
// traditional sum of squares and products method of calculation
// with subtraction of means - less subject to round-off error
// than test1
// make matrix of predictor values
Matrix X(nobs,npred);
// load x1 and x2 into X
// [use << rather than = when loading arrays]
X.Column(1) << x1; X.Column(2) << x2;
// vector of Y values
ColumnVector Y(nobs); Y << y;
// make vector of 1s
ColumnVector Ones(nobs); Ones = 1.0;
// calculate means (averages) of x1 and x2 [ .t() takes transpose]
RowVector M = Ones.t() * X / nobs;
// and subtract means from x1 and x1
Matrix XC(nobs,npred);
XC = X - Ones * M;
// do the same to Y [use Sum to get sum of elements]
ColumnVector YC(nobs);
Real m = Sum(Y) / nobs; YC = Y - Ones * m;
// form sum of squares and product matrix
// [use << rather than = for copying Matrix into SymmetricMatrix]
SymmetricMatrix SSQ; SSQ << XC.t() * XC;
// calculate estimate
// [bracket last two terms to force this multiplication first]
// [ .i() means inverse, but inverse is not explicity calculated]
ColumnVector A = SSQ.i() * (XC.t() * YC);
// calculate estimate of constant term
// [AsScalar converts 1x1 matrix to Real]
Real a = m - (M * A).AsScalar();
// Get variances of estimates from diagonal elements of inverse of SSQ
// [ we are taking inverse of SSQ - we need it for finding D ]
Matrix ISSQ = SSQ.i(); DiagonalMatrix D; D << ISSQ;
ColumnVector V = D.AsColumn();
Real v = 1.0/nobs + (M * ISSQ * M.t()).AsScalar();
// for calc variance of const
// Calculate fitted values and residuals
int npred1 = npred+1;
ColumnVector Fitted = X * A + a;
ColumnVector Residual = Y - Fitted;
Real ResVar = Residual.SumSquare() / (nobs-npred1);
// Get diagonals of Hat matrix (an expensive way of doing this)
Matrix X1(nobs,npred1); X1.Column(1)<<Ones; X1.Columns(2,npred1)<<X;
DiagonalMatrix Hat; Hat << X1 * (X1.t() * X1).i() * X1.t();
// print out answers
cout << "\nEstimates and their standard errors\n\n";
cout.setf(ios::fixed, ios::floatfield);
cout << setw(11) << setprecision(5) << a << " ";
cout << setw(11) << setprecision(5) << sqrt(v*ResVar) << endl;
// make vector of standard errors
ColumnVector SE(npred);
for (int i=1; i<=npred; i++) SE(i) = sqrt(V(i)*ResVar);
// use concatenation function to form matrix and use matrix print
// to get two columns
cout << setw(11) << setprecision(5) << (A | SE) << endl;
cout << "\nObservations, fitted value, residual value, hat value\n";
cout << setw(9) << setprecision(3) <<
(X | Y | Fitted | Residual | Hat.AsColumn());
cout << "\n\n";
}
示例4: trymatd
void trymatd()
{
Tracer et("Thirteenth test of Matrix package");
Tracer::PrintTrace();
Matrix X(5,20);
int i,j;
for (j=1;j<=20;j++) X(1,j) = j+1;
for (i=2;i<=5;i++) for (j=1;j<=20; j++) X(i,j) = (long)X(i-1,j) * j % 1001;
SymmetricMatrix S; S << X * X.t();
Matrix SM = X * X.t() - S;
Print(SM);
LowerTriangularMatrix L = Cholesky(S);
Matrix Diff = L*L.t()-S; Clean(Diff, 0.000000001);
Print(Diff);
{
Tracer et1("Stage 1");
LowerTriangularMatrix L1(5);
Matrix Xt = X.t(); Matrix Xt2 = Xt;
QRZT(X,L1);
Diff = L - L1; Clean(Diff,0.000000001); Print(Diff);
UpperTriangularMatrix Ut(5);
QRZ(Xt,Ut);
Diff = L - Ut.t(); Clean(Diff,0.000000001); Print(Diff);
Matrix Y(3,20);
for (j=1;j<=20;j++) Y(1,j) = 22-j;
for (i=2;i<=3;i++) for (j=1;j<=20; j++)
Y(i,j) = (long)Y(i-1,j) * j % 101;
Matrix Yt = Y.t(); Matrix M,Mt; Matrix Y2=Y;
QRZT(X,Y,M); QRZ(Xt,Yt,Mt);
Diff = Xt - X.t(); Clean(Diff,0.000000001); Print(Diff);
Diff = Yt - Y.t(); Clean(Diff,0.000000001); Print(Diff);
Diff = Mt - M.t(); Clean(Diff,0.000000001); Print(Diff);
Diff = Y2 * Xt2 * S.i() - M * L.i();
Clean(Diff,0.000000001); Print(Diff);
}
ColumnVector C1(5);
{
Tracer et1("Stage 2");
X.ReSize(5,5);
for (j=1;j<=5;j++) X(1,j) = j+1;
for (i=2;i<=5;i++) for (j=1;j<=5; j++)
X(i,j) = (long)X(i-1,j) * j % 1001;
for (i=1;i<=5;i++) C1(i) = i*i;
CroutMatrix A = X;
ColumnVector C2 = A.i() * C1; C1 = X.i() * C1;
X = C1 - C2; Clean(X,0.000000001); Print(X);
}
{
Tracer et1("Stage 3");
X.ReSize(7,7);
for (j=1;j<=7;j++) X(1,j) = j+1;
for (i=2;i<=7;i++) for (j=1;j<=7; j++)
X(i,j) = (long)X(i-1,j) * j % 1001;
C1.ReSize(7);
for (i=1;i<=7;i++) C1(i) = i*i;
RowVector R1 = C1.t();
Diff = R1 * X.i() - ( X.t().i() * R1.t() ).t(); Clean(Diff,0.000000001);
Print(Diff);
}
{
Tracer et1("Stage 4");
X.ReSize(5,5);
for (j=1;j<=5;j++) X(1,j) = j+1;
for (i=2;i<=5;i++) for (j=1;j<=5; j++)
X(i,j) = (long)X(i-1,j) * j % 1001;
C1.ReSize(5);
for (i=1;i<=5;i++) C1(i) = i*i;
CroutMatrix A1 = X*X;
ColumnVector C2 = A1.i() * C1; C1 = X.i() * C1; C1 = X.i() * C1;
X = C1 - C2; Clean(X,0.000000001); Print(X);
}
{
Tracer et1("Stage 5");
int n = 40;
SymmetricBandMatrix B(n,2); B = 0.0;
for (i=1; i<=n; i++)
{
B(i,i) = 6;
if (i<=n-1) B(i,i+1) = -4;
if (i<=n-2) B(i,i+2) = 1;
}
B(1,1) = 5; B(n,n) = 5;
SymmetricMatrix A = B;
ColumnVector X(n);
X(1) = 429;
for (i=2;i<=n;i++) X(i) = (long)X(i-1) * 31 % 1001;
X = X / 100000L;
// the matrix B is rather ill-conditioned so the difficulty is getting
// good agreement (we have chosen X very small) may not be surprising;
// maximum element size in B.i() is around 1400
ColumnVector Y1 = A.i() * X;
LowerTriangularMatrix C1 = Cholesky(A);
ColumnVector Y2 = C1.t().i() * (C1.i() * X) - Y1;
Clean(Y2, 0.000000001); Print(Y2);
UpperTriangularMatrix CU = C1.t().i();
//.........这里部分代码省略.........