/* * * Copyright 1996 Numerical Algorithms Group * * Include file for NAG C Library g03 Chapter * * Mark 5, 1997. * Mark 6 revised. IER-3042 (July 2000). */ #ifndef NAGG03 #define NAGG03 //#importdll "ONAG" // NAG DLL prepared by OriginLab #pragma dll(ONAG) #include /* begin proto */ /** g03aac performs a principal component analysis on a data matrix; both the principal component loadings and the principal component scores are returned. Example 1: Assume we have four Matrix, "Matrix1"'s size is 10X3 and stores 30 number; "Matrix2"'s size is 3X6; "Matrix3"'s size is 3X3, and "Matrix4"'s size is 12X3. We will store the value to "Matrix2", "Matrix3" and "Matrix4". matrix p; p.SetSize(3,3); double s[3]; matrix e; e.SetSize(3,6); matrix v; v.SetSize(12,3); double wt[12]; double *wtptr; wtptr = NULL; int isx[3] = {1, 1, 1}; int nvar = 3; int tdx = 3; int tde = 6; int tdp = 3; int tdv = 3; int i, j, m = 3, n = 10; Nag_PrinCompMat pcmatrix; Nag_PrinCompScores scores; Matrix xx("Matrix1"); matrix x = xx; char weight, mat, std; mat = 'V'; std = 'E'; weight = 'U'; if (mat == 'C') pcmatrix = Nag_MatCorrelation; else if (mat == 'S') pcmatrix = Nag_MatStandardised; else if (mat == 'U') pcmatrix = Nag_MatSumSq; else pcmatrix = Nag_MatVarCovar; if (std == 'S') scores = Nag_ScoresStand; else if (std == 'U') scores = Nag_ScoresNotStand; else if (std == 'Z') scores = Nag_ScoresUnitVar; else scores = Nag_ScoresEigenval; nag_mv_prin_comp(pcmatrix, scores, n, m, x, tdx, isx, s, wtptr, nvar, e, tde, p, tdp, v, tdv); //put the result to the matrix. //first column is Eigenvalues, second column is Percentage variation, //third column is Cumulative variation, fourth column is Chisp, //fifth column is DF and sixth column is Sig. Matrix ee("Matrix2"); ee = e; Matrix pp("Matrix3"); //Eigenvalues pp = p; Matrix vv("Matrix4"); //Principal component scores vv = v; Example 2: A data set is taken from Cooley and Lohnes (1971), it consists of ten observations on three variables. The unweighted principal components based on the variance-covariance matrix are computed and unstandardised principal component scores requested. void test_nag_mv_prin_comp() { matrix p; p.SetSize(3,3); double s[3]; matrix e; e.SetSize(3,6); matrix v; v.SetSize(12,3); matrix x = {{7.0, 4.0, 3.0}, {4.0, 1.0, 8.0}, {6.0, 3.0, 5.0}, {8.0, 6.0, 1.0}, {8.0, 5.0, 7.0}, {7.0, 2.0, 9.0}, {5.0, 3.0, 3.0}, {9.0, 5.0, 8.0}, {7.0, 4.0, 5.0}, {8.0, 2.0, 2.0}}; double wt[12]; double *wtptr; int isx[3]; int nvar = 3; int tdx = 3; int tde = 6; int tdp = 3; int tdv = 3; int i, j, m = 3, n = 10; Nag_PrinCompMat pcmatrix; Nag_PrinCompScores scores; char weight, mat, std; for(i = 0; i < 3; i++) isx[i] = 1; mat = 'V'; std = 'E'; weight = 'U'; if (mat == 'C') pcmatrix = Nag_MatCorrelation; else if (mat == 'S') pcmatrix = Nag_MatStandardised; else if (mat == 'U') pcmatrix = Nag_MatSumSq; else pcmatrix = Nag_MatVarCovar; if (std == 'S') scores = Nag_ScoresStand; else if (std == 'U') scores = Nag_ScoresNotStand; else if (std == 'Z') scores = Nag_ScoresUnitVar; else scores = Nag_ScoresEigenval; wtptr = NULL; nag_mv_prin_comp(pcmatrix, scores, n, m, x, tdx, isx, s, wtptr, nvar, e, tde, p, tdp, v, tdv); printf("Eigenvalues Percentage Cumulative Chisq DF Sig\n"); printf(" variation variation\n\n"); for (i = 0; i < nvar; ++i) { for (j = 0; j < 6; ++j) printf("%11.4f",e[i][j]); printf("\n"); } printf("\nEigenvalues\n\n"); for (i = 0; i < nvar; ++i) { for (j = 0; j < nvar; ++j) printf("%9.4f",p[i][j]); printf("\n"); } printf("\nPrincipal component scores \n\n"); for (i = 0; i < n; ++i) { printf("%2ld", i+1); for (j = 0; j < nvar; ++j) printf("%9.3f", v[i][j]); printf("\n"); } } The output is following: Eigenvalues Percentage Cumulative Chisq DF Sig variation variation 8.2739 0.6515 0.6515 8.6127 5.0000 0.1255 3.6761 0.2895 0.9410 4.1183 2.0000 0.1276 0.7499 0.0590 1.0000 0.0000 0.0000 0.0000 Eigenvalues 0.1376 0.6990 0.7017 0.2505 0.6609 -0.7075 -0.9583 0.2731 -0.0842 Principal component scores 1 2.151 -0.173 -0.107 2 -3.804 -2.887 -0.510 3 -0.153 -0.987 -0.269 4 4.707 1.302 -0.652 5 -1.294 2.279 -0.449 6 -4.099 0.144 0.803 7 1.626 -2.232 -0.803 8 -2.114 3.251 0.168 9 0.235 0.373 -0.275 10 2.746 -1.069 2.094 Return: The function returns NAG error code or 0 if no error. 70: On entry, parameter pcmatrix had an illegal value. On entry, parameter scores had an illegal value. 11: On entry, m must not be less than 1: m = _value_. On entry, n must not be less than 2: n = _value_. On entry, nvar must not be less than 1: nvar = _value_. On entry, tde must not be less than 6: tde = _value_. 17: On entry, tdx = _value_ while m = _value_. These parameters must satisfy tdx = m. On entry, tdv = _value_ while nvar = _value_. These parameters must satisfy tdv = nvar. On entry, tdp = _value_ while nvar = _value_. These parameters must satisfy tdp = nvar. 19: On entry, nvar = _value_ while m = _value_. These parameters must satisfy nvar = m. 20: On entry, nvar = _value_ while n = _value_. These parameters must satisfy nvar < n. 500: On entry, wt[_value_] = _value_. Constraint: when referenced, all elements of wt must be non-negative. 501: The number of variables, nvar in the analysis = _value_, while the number of variables included in the analysis via array isx = _value_. Constraint: these two numbers must be the same. 504: On entry, the standardisation element s[_value_] = _value_, while the variable to be included isx[_value_] = _value_. Constraint: when a variable is to included, the standardisation element must be positive. 503: With weighted data, the e.ective number of observations given by the sum of weights = _value_, while the number of variables included in the analysis, nvar = _value_. Constraint: effective number of observations > nvar + 1. 427: The singular value decomposition has failed to converge. This is an unlikely error exit. 505: All eigenvalues/singular values are zero. This will be caused by all the variables being constant. 73: Memory allocation failed. 74: An internal error has occurred in this function. Check the function call and array sizes. If the function call is correct, please consult NAG for assistance. successfully call of the nag_mv_prin_comp function. */ int nag_mv_prin_comp( Nag_PrinCompMat pcmatrix, Nag_PrinCompScores scores, int n, // the number of observations int m, // the number of variables in the data matrix double x[], // contain the ith observation for the jth variable, for i = 1, 2, . . ., n; j = 1, 2, . . .,m. int tdx, // the last dimension of the array x. int isx[], // indicates whether or not the jth variable is to be included in the analysis. double s[], // the standardizations to be used, if any. double w[], // contain the weights to be used in the principal component analysis. int nvar, // the number of variables in the principal component analysis. double e[], // the statistics of the principal component analysis. int tde, // the last dimension of the array e. double p[], // the first nvar columns of p contain the principal component loadings. The jth column of p contains the nvar coefficients for the jth principal component. int tdp, // the last dimension of the array p double v[], // the first nvar columns of v contain the principal component scores. The jth column of v contains the n scores for the jth principal component. int tdv // the last dimension of the array v. ); // Eigenvalues, Percentage Variation, Cumulative Variation, Principal Component Scores...... /** g03acc performs a canonical variate (canonical descrimination) analysis. Example: A sample of nine observations, each consisting of three variables plus group indicator, is read in. There are three groups. An unweighted canonical variate analysis is performed. void test_nag_mv_canon_var() { matrix e; e.SetSize(3,6); matrix x = {{13.3, 10.6, 21.2}, {13.6, 10.2, 21.0}, {14.2, 10.7, 21.1}, {13.4, 9.4, 21.0}, {13.2, 9.6, 20.1}, {13.9, 10.4, 19.8}, {12.9, 10.0, 20.5}, {12.2, 9.9, 20.7}, {13.9, 11.0, 19.1}}; double wt[9]; matrix cvm; cvm.SetSize(3,3); double tol; matrix cvx; cvx.SetSize(3,3); int i, j, m = 3, n = 9; int ng = 3; int nx = 3; int nig[3], ncv; int ing[9] = {1, 2, 3, 1, 2, 3, 1, 2, 3}; int irx, isx[6]; int tdx = 3; int tdc = 3; int tde = 6; char wtchar = 'U'; Nag_Weightstype weight; weight = Nag_NoWeights; for (j = 0; j < m; ++j) isx[j] = 1; tol = 1e-6; nag_mv_canon_var(weight, n, m, x, tdx, isx, nx, ing, ng, wt, nig, cvm, tdc, e, tde, &ncv, cvx, tdc, tol, &irx); printf("%s%2ld\n\n","Rank of x = ",irx); printf(" Canonical Eigenvalues Percentage ChiSQ DF SIG\n"); printf(" Correlations Variation\n"); for (i = 0; i < ncv; ++i) { for (j = 0; j < 6; ++j) printf("%12.4f",e[i][j]); printf("\n"); } printf("\nCanonical Coefficients for X\n"); for (i = 0; i < nx; ++i) { for (j = 0; j < ncv; ++j) printf("%9.4f",cvx[i][j]); printf("\n"); } printf("\nCanonical variate means\n"); for (i = 0; i < ng; ++i) { for (j = 0; j < ncv; ++j) printf("%9.4f",cvm[i][j]); printf("\n"); } } The output is following: Rank of x = 3 Canonical Eigenvalues Percentage CHISQ DF SIG Correlations Variation 0.8826 3.5238 0.9795 7.9032 6.0000 0.2453 0.2623 0.0739 0.0205 0.3564 2.0000 0.8368 Canonical Coefficients for X -1.7070 0.7277 -1.3481 0.3138 0.9327 1.2199 Canonical variate means 0.9841 0.2797 1.1805 -0.2632 -2.1646 -0.0164 Return: The function returns NAG error code or 0 if no error. 70: On entry, parameter weight had an illegal value. 11: On entry, nx must not be less than 1: nx = _value_. On entry, ng must not be less than 2: ng = _value_. On entry, tde must not be less than 6: tde = _value_. 5: On entry, tol must not be less than 0.0: tol = _value_. 17: On entry, m = _value_ while nx = _value_. These parameters must satisfy m = nx. On entry, tdx = _value_ while m = _value_. These parameters must satisfy tdx = m. On entry, tdcvx = _value_ while ng = _value_. These parameters must satisfy tdcvx = ng-1. On entry, tdcvm = _value_ while nx = _value_. These parameters must satisfy tdcvm = nx. 29: On entry, n = _value_, nx = _value_ and ng = _value_. These parameters must satisfy n = nx + ng. 112: On entry, ing[_value_] = _value_, ng = _value_. Constraint: 1 = ing[i - 1] = ng, i = 1, 2, . . ., n. 510: The wt array argument must not be NULL when the weight argument indicates weights. 500: On entry, wt[_value_] = _value_. Constraint: When referenced, all elements of wt must be non-negative. 501: The number of variables, nx in the analysis = _value_, while number of variables included in the analysis via array isx = _value_. Constraint: these two numbers must be the same. 427: The singular value decomposition has failed to converge. This is an unlikely error exit. 506: A canonical correlation is equal to one. This will happen if the variables provide an exact indication as to which group every observation is allocated. 507: Either the e.ective number of groups is less than two or the e.ective number of groups plus the number of variables, nx is greater than the the e.ective number of observations. 508: The rank of the variables is zero. This will happen if all the variables are constants. 73: Memory allocation failed. 74: An internal error has occurred in this function. Check the function call and array sizes. If the function call is correct, please consult NAG for assistance. successfully call of the nag_mv_canon_var function. */ int nag_mv_canon_var( Nag_Weightstype weights, //indicates the type of weights to be used in the analysis. int n, // the number of observations. int m, // the total number of variables. double x[], // contain the ith observation for the jth variable, for i = 1, 2, . . ., n; j = 1, 2, . . .,m. int tdx, // the last dimension of the array x. int isx[], // indicates whether or not the jth variable is to be included in the analysis. int nx, // the number of variables in the analysis. int ing[], // indicates which group the ith observation is in, for i = 1, 2, . . ., n. int ng, // The number of groups. double *wtptr, // if weight = Nag Weightsfreq or Nag Weightsvar then the elements of wt must contain the weights to be used in the analysis. int nig[], // the number of observations in group j, for j = 1, 2, . . . , ng. double cvm[], // the mean of the jth canonical variate for the ith group, for i = 1, 2, . . . , ng; j = 1, 2, . . . , l; int tdcvm, // the last dimension of the array cvm. double e[], // the statistics of the canonical variate analysis. int tde, // the last dimension of the array e. int *ncv, // the number of canonical variates. double cvx[], // the canonical variate loadings. int tdcvx, // the last dimension of the array cvx. double tol, // the value of tol is used to decide if the variables are of full rank. int *irankx // the rank of the dependent variables. ); // Canonical Coeffiecients, Canonical Correlation, Eigenvalues, Canonical Variate Means.... /** g03adc performs canonical correlation analysis upon input data matrices. Example: A sample of nine observations with two variables in each set is read in. The second and third variables are x variables while the first and last are y variables. Canonical variate analysis is performed. void test_nag_mv_canon_corr() { matrix e, cvx, cvy; e.SetSize(2,6); cvx.SetSize(2,2); cvy.SetSize(2,2); matrix z = {{80.0, 58.4, 14.0, 21.0}, {75.0, 59.2, 15.0, 27.0}, {78.0, 60.3, 15.0, 27.0}, {75.0, 57.4, 13.0, 22.0}, {79.0, 59.5, 14.0, 26.0}, {78.0, 58.1, 14.5, 26.0}, {75.0, 58.0, 12.5, 23.0}, {64.0, 55.5, 11.0, 22.0}, {80.0, 59.2, 12.5, 22.0}}; double wt[9]; double tol; double *wtptr; int i, j, m = 4, n = 9; int ix = 2, iy = 2; int nx, ny; int ncv; int isz[4] = {-1, 1, 1, -1}; int tdz= 4; int tde=6; int tdc=2; char weight = 'U'; wtptr = NULL; tol = 1e-6; nx = ix; ny= iy; nag_mv_canon_corr(n, m, z, tdz, isz, nx, ny, wtptr, e, tde, &ncv, cvx, tdc, cvy, tdc, tol); printf("\n%s%2ld%s%2ld\n\n","Rank of x = ",nx," Rank of y= ",ny); printf(" Canonical Eigenvalues Percentage Chisq DF Sig\n"); printf(" correlations variation\n"); for (i = 0; i < ncv; ++i) { for (j = 0; j < 6; ++j) printf("%12.4f",e[i][j]); printf("\n"); } printf("\nCanonical coefficients for x\n"); for (i = 0; i < ix; ++i) { for (j = 0; j < ncv; ++j) printf("%9.4f",cvx[i][j]); printf("\n"); } printf("\nCanonical coefficients for y\n"); for (i = 0; i < iy ;++i) { for (j = 0; j < ncv; ++j) printf("%9.4f",cvy[i][j]); printf("\n"); } } The output is following: Rank of x = 2 Rank of y= 2 Canonical Eigenvalues Percentage Chisq DF Sig correlations variation 0.9570 10.8916 0.9863 14.3914 4.0000 0.0061 0.3624 0.1512 0.0137 0.7744 1.0000 0.3789 Canonical coefficients for x -0.4261 1.0337 -0.3444 -1.1136 Canonical coefficients for y -0.1415 0.1504 -0.2384 -0.3424 Return: The function returns NAG error code or 0 if no error. 11: On entry, nx must not be less than 1: nx = _value_. On entry, ny must not be less than 1: ny = _value_. On entry, tde must not be less than 6: tde = _value_. 5: On entry, tol must not be less than 0.0: tol = _value_. 17: On entry, tdz = _value_ while m = _value_. These parameters must satisfy tdz = m. 29: On entry, m = _value_, nx = _value_ and ny = _value_. These parameters must satisfy m = nx + ny. On entry, n = _value_, nx = _value_ and ny = _value_. These parameters must satisfy n > nx + ny. On entry, tdcvx = _value_, nx = _value_ and ny = _value_. These parameters must satisfy tdcvx = min(nx,ny). On entry, tdcvy = _value_, nx = _value_ and ny = _value_. These parameters must satisfy tdcvy = min(nx,ny). 500: On entry, wt[_value_] = _value_. Constraint: When referenced, all elements of wt must be non-negative. 501: The number of variables, nx in the analysis = _value_, while the number of x variables included in the analysis via array isz = _value_. Constraint: these two numbers must be the same. The number of variables, ny in the analysis = _value_, while the number of y variables included in the analysis via array isz = _value_. Constraint: these two numbers must be the same. 503: With weighted data, the e.ective number of observations given by the sum of weights = _value_, while number of variables included in the analysis, nx + ny = _value_. Constraint: E.ective number of observations = nx + ny + 1. 427: The singular value decomposition has failed to converge. This is an unlikely error exit. 506: A canonical correlation is equal to one. This will happen if the x and y variables are perfectly correlated. 509: The rank of the x matrix or the rank of the y matrix is zero. This will happen if all the x and y variables are constants. 73: Memory allocation failed. 74: An internal error has occurred in this function. Check the function call and array sizes. If the function call is correct, please consult NAG for assistance. successfully call of the nag_mv_canon_corr function. */ int nag_mv_canon_corr( int n, // the number of observations. int m, // the total number of variables. double z[], // contain the ith observation for the jth variable, for i = 1, 2, . . ., n; j = 1, 2, . . .,m. int tdz, // the last dimension of the array z. int isz[], // indicates whether or not the jth variable is to be included in the analysis and to which set of variables it belongs. int nx, // the number of x variables in the analysis. int ny, // the number of y variables in the analysis. double wt[], // the weights to be used in the analysis. double e[], // the statistics of the canonical variate analysis. int tde, // the last dimension of the array e. int *ncv, // the number of canonical correlations. double cvx[], // the canonical variate loadings for the x variables. int tdcvx, // the last dimension of the array cvx. double cvy[], // the canonical variate loadings for the y variables. int tdcvy, // the last dimension of the array cvy. double tol ); // Canonical coefficients, Eigenvalues, Percentage variation, Canonical coefficients... /** g03bac computes orthogonal rotations for a matrix of loadings using a generalized orthomax criterion. Example: The example is taken from page 75 of Lawley and Maxwell (1971). The results from a factor analysis of ten variables using three factors are input and rotated using varimax rotations without standardising rows. void test_nag_mv_orthomax() { double g = 1.0; matrix r, flr; r.SetSize(3,3); flr.SetSize(10,3); matrix fl = {{0.788, -0.152, -0.352}, {0.874, 0.381, 0.041}, {0.814, -0.043, -0.213}, {0.798, -0.170, -0.204}, {0.641, 0.070, -0.042}, {0.755, -0.298, 0.067}, {0.782, -0.221, 0.028}, {0.767, -0.091, 0.358}, {0.733, -0.384, 0.229}, {0.771, -0.101, 0.071}}; double acc = 0.00001; int iter, nvar; int i, j, k; int maxit = 20; int tdf = 3; int tdr = 3; char char_stand; Nag_RotationLoading stand; nvar = 10; k = 3; char_stand = 'U'; if (char_stand == 'S') stand = Nag_RoLoadStand; else stand = Nag_RoLoadNotStand; nag_mv_orthomax(stand, g, nvar, k, fl, tdf, flr, r, tdr, acc, maxit, &iter); printf("\n Rotated factor loadings\n\n"); for (i = 0; i < nvar; ++i) { for (j = 0; j < k ;++j) printf(" %8.3f",flr[i][j]); printf("\n"); } printf("\n Rotation matrix\n\n"); for (i = 0; i < k ;++i) { for (j = 0; j < k ;++j) printf(" %8.3f",r[i][j]); printf("\n"); } } The output is following: Rotated factor loadings 0.329 -0.289 -0.759 0.849 -0.273 -0.340 0.450 -0.327 -0.633 0.345 -0.397 -0.657 0.453 -0.276 -0.370 0.263 -0.615 -0.464 0.332 -0.561 -0.485 0.472 -0.684 -0.183 0.209 -0.754 -0.354 0.423 -0.514 -0.409 Rotation matrix 0.633 -0.534 -0.560 0.758 0.573 0.311 0.155 -0.622 0.768 Return: The function returns NAG error code or 0 if no error. 70: On entry, parameter stand had an illegal value. 11: On entry, k must not be less than 2: k = _value_. 12: On entry, maxit must not be less than or equal to 0 : maxit = _value_. 5: On entry, g must not be less than 0.0: g = _value_. On entry, acc must not be less than 0.0: acc = _value_. 17: On entry, nvar = _value_ while k = _value_. These parameters must satisfy nvar = k. On entry, tdf = _value_ while k = _value_. These parameters must satisfy tdf = k. On entry, tdr = _value_ while k = _value_. These parameters must satisfy tdr = k. 427: The singular value decomposition has failed to converge. This is an unlikely error exit. 511: The algorithm to .nd R has failed to reach the required accuracy in the given number of iterations, _value_. Try increasing acc or increasing maxit. The returned solution should be a reasonable approximation. 73: Memory allocation failed. 74: An internal error has occurred in this function. Check the function call and array sizes. If the function call is correct, please consult NAG for assistance. successfully call of the nag_mv_orthomax function. */ int nag_mv_orthomax( Nag_RotationLoading stand, double g, // the criterion constant. int nvar, // The number of original variables. int k, // The number of derived variates or factors. double fl[], // the matrix of loadings. int tdf, // the last dimension of the arrays fl and flr. double flr[], // the rotated matrix of loadings. double r[], // the matrix of rotations. int tdr, // the last dimension of the array r. double acc, // indicates the accuracy required. int maxit, // the maximum number of iterations. int *iter // the number of iterations performed. ); // Rotation factor loading, Rotation matrix. /** g03bcc computes Procrustes rotations in which an orthogonal rotation is found so that a transformed matrix best matches a target matrix. Example: Three points representing the vertices of a triangle in two dimensions are input. The points are translated and rotated to match the triangle given by (0,0),(1,0),(0,2) and scaling is applied after rotation. void test_nag_mv_procustes() { matrix r, yhat; r.SetSize(2,2); yhat.SetSize(3,2); matrix x = {{0.63, 0.58}, {1.36, 0.39}, {1.01, 1.76}}; matrix y = {{0.0, 0.0}, {1.0, 0.0}, {0.0, 2.0}}; double res[3]; double alpha; double rss; int i, j, m = 2, n = 3; int tdx = 2; int tdr = 2; int tdy = 2; char char_scale; char char_stand; Nag_TransNorm stand; Nag_RotationScale scale; char_stand = 'C'; char_scale = 'S'; if (char_stand == 'N') stand = Nag_NoTransNorm; else if (char_stand == 'Z') stand = Nag_Orig; else if (char_stand == 'C') stand = Nag_OrigCentroid; else if (char_stand == 'U') stand = Nag_Norm; else if (char_stand == 'S') stand = Nag_OrigNorm; else if (char_stand == 'M') stand = Nag_OrigNormCentroid; if (char_scale == 'S') scale = Nag_LsqScale; else if (char_scale == 'U') scale = Nag_NotLsqScale; nag_mv_procustes(stand, scale, n, m, x, tdx, y, tdy, yhat, r, tdr, &alpha, &rss, res); printf("\n Rotation matrix\n\n"); for (i = 0; i < m; ++i) { for (j = 0; j < m; ++j) printf(" %7.3f ",r[i][j]); printf("\n"); } if (char_scale == 'S' || char_scale == 's') printf("\n%s%10.3f\n"," Scale factor = ",alpha); printf("\n Target matrix \n\n"); for (i = 0; i < n; ++i) { for (j = 0; j < m; ++j) printf(" %7.3f ",y[i][j]); printf("\n"); } printf("\n Fitted matrix\n\n"); for (i = 0; i < n; ++i) { for (j = 0; j < m; ++j) printf(" %7.3f ",yhat[i][j]); printf("\n"); } printf("\n%s%10.3f\n","RSS = ",rss); } The output is following: Rotation matrix 0.967 0.254 -0.254 0.967 Scale factor = 1.556 Target matrix 0.000 0.000 1.000 0.000 0.000 2.000 Fitted matrix -0.093 0.024 1.080 0.026 0.013 1.950 RSS = 0.019 Return: The function returns NAG error code or 0 if no error. 70: On entry, parameter stand had an illegal value. On entry, parameter pscale had an illegal value. 11: On entry, n must not be less than 1: n = _value_. On entry, m must not be less than 1: m = _value_. 17: On entry, tdx = _value_ while m = _value_. These parameters must satisfy tdx = m. On entry, tdy = _value_ while m = _value_. These parameters must satisfy tdy = m. On entry, tdr = _value_ while m = _value_. These parameters must satisfy tdr = m. 19: On entry, m = _value_ while n = _value_. These parameters must satisfy m = n. 427: The singular value decomposition has failed to converge. This is an unlikely error exit. 512: On entry, either x or y contains only zero-points (possibly after translation) when normalization is to be applied. 513: The fitted matrix Y-hat, contains only zero-points when least-squares scaling is applied. 73: Memory allocation failed. 74: An internal error has occurred in this function. Check the function call and array sizes. If the function call is correct, please consult NAG for assistance. successfully call of the nag_mv_orthomax function. */ int nag_mv_procustes( Nag_TransNorm stand, Nag_RotationScale pscale, int n, // the number of points. int m, // the number of dimensions. double x[], // the matrix to be rotated. int tdx, // the last dimension of the array x. double y[], // the target matrix. int tdy, // the last dimension of the arrays y. double yhat[], // the fitted matrix double r[], // the matrix of rotations. int tdr, // the last dimension of the array r. double *alpha, // if pscale = Nag LsqScale the scaling factor, a; otherwise alpha is not set double *rss, // the residual sum of squares. double res[] // the residuals. ); // Rotation matrix, Target matrix, Fitted matrix, RSS, Scale factor. /** g03cac computes the maximum likelihood estimates of the parameters of a factor analysis model. Either the data matrix or a correlation/covariance matrix may be input. Factor loadings, communalities and residual correlations are returned. Example: The example is taken from Lawley and Maxwell (1971). The correlation matrix for nine variables is input and the parameters of a factor analysis model with three factors are estimated and printed. // Note: This function needs the support of struct. It has not been done yet! void test_nag_mv_factor() { matrix fl, x; fl.SetSize(9,9); x.SetSize(9,9); double com[9], e[9]; double psi[9], res[36], stat[4], wt[9]; double *wtptr=NULL; double eps; int nfac, nvar; int i, j, l, m, n; int isx[9]; int tdfl = 9, tdx = 9; char weight, char_mat; Nag_FacMat mat; Nag_E04_Opt options; Vprintf("g03cac Example Program Results\n\n"); Vscanf("%*[^\n]"); Vscanf("%s",char_mat); Vscanf("%s",weight); Vscanf("%ld",&n); Vscanf("%ld",&m); Vscanf("%ld",&nvar); Vscanf("%ld",&nfac); if (m <= MMAX && (*char_mat == 'C' || n <= NMAX )) { if (*char_mat == 'C') { for (i = 0; i < m; ++i) { for (j = 0; j < m; ++j) Vscanf("%lf",&x[i][j]); } } else { if (*weight == 'W' || *weight == 'w') { for (i = 0; i < n; ++i) { for (j = 0; j < m; ++j) Vscanf("%lf",&x[i][j]); Vscanf("%lf",&wt[i]); } wtptr = wt; } else { for (i = 0; i < n; ++i) { for (j = 0; j < m; ++j) Vscanf("%lf",&x[i][j]); } } } for (j = 0; j < m; ++j) Vscanf("%ld",&isx[j]); 3.g03cac.6 [NP3275/5/pdf] g03 - Multivariate Methods g03cac if (*char_mat == 'D') { mat = Nag_DataCorr; } else if (*char_mat == 'S') { mat = Nag_DataCovar; } else if (*char_mat == 'C') { mat = Nag_MatCorr_Covar; } e04xxc(&options); options.max_iter = 500; options.optim_tol = 1e-2; eps = 1e-5; g03cac(mat, n, m, (double *)x, tdx, nvar, isx, nfac, wtptr, e, stat, com, psi, res, (double *)fl, tdfl, &options, eps, NAGERR_DEFAULT); Vprintf("\nEigenvalues\n\n"); for (j = 0; j < m; ++j) { Vprintf(" %12.4e%s",e[j], (j+1)%6==0 ? "\n" : ""); } Vprintf("\n\n%s%6.3f\n"," Test Statistic = ",stat[1]); Vprintf("%s%6.3f\n"," df = ",stat[2]); Vprintf("%s%6.3f\n\n","Significance level = ",stat[3]); Vprintf("Residuals\n\n"); l = 1; for (i = 1; i <= nvar-1; ++i) { for (j = l; j <= l+i-1; ++j) Vprintf(" %8.3f",res[j-1]); Vprintf("\n"); l += i; } Vprintf("\nLoadings, Communalities and PSI\n\n"); for (i = 0; i < nvar; ++i) { for (j = 0; j < nfac; ++j) Vprintf(" %8.3f",fl[i][j]); Vprintf("%8.3f%8.3f\n", com[i], psi[i]); } exit(EXIT_SUCCESS); } else { Vprintf("Incorrect input value of n or m.\n"); exit(EXIT_FAILURE); } } The output is following: Parameters to e04lbc -------------------- Number of variables........... 9 optim_tol............... 1.00e-02 linesearch_tol.......... 9.00e-01 step_max................ 2.70e+01 max_iter................ 500 print_level.........Nag_Soln_Iter machine precision....... 1.11e-16 deriv_check............. FALSE outfile................. stdout Memory allocation: state................... User hesl.................... User hesd................... User Iterations performed = 0, function evaluations = 1 Criterion = 8.635756e-02 Variable Standardized Communalities 1 0.5755 2 0.5863 3 0.4344 4 0.7496 5 0.6203 6 0.7329 7 0.6061 8 0.4053 9 0.7104 Iterations performed = 1, function evaluations = 3 Criterion = 3.603203e-02 Variable Standardized Communalities 1 0.5517 2 0.5800 3 0.3936 4 0.7926 5 0.6140 6 0.8254 7 0.6052 8 0.5076 9 0.7569 Iterations performed = 2, function evaluations = 4 Criterion = 3.502097e-02 Variable Standardized Communalities 1 0.5496 2 0.5731 3 0.3838 4 0.7875 5 0.6200 6 0.8238 7 0.6006 8 0.5349 9 0.7697 Iterations performed = 3, function evaluations = 5 Criterion = 3.501729e-02 Variable Standardized Communalities 1 0.5495 2 0.5729 3 0.3835 4 0.7877 5 0.6195 6 0.8231 7 0.6005 8 0.5384 9 0.7691 Eigenvalues 1.5968e+01 4.3577e+00 1.8474e+00 1.1560e+00 1.1190e+00 1.0271e+00 9.2574e-01 8.9508e-01 8.7710e-01 Test Statistic = 7.149 df = 12.000 Significance level = 0.848 Residuals 0.000 -0.013 0.022 0.011 -0.005 0.023 -0.010 -0.019 -0.016 0.003 -0.005 0.011 -0.012 -0.001 -0.001 0.015 -0.022 -0.011 0.002 0.029 -0.012 -0.001 -0.011 0.013 0.005 -0.006 -0.001 0.003 -0.006 0.010 -0.005 -0.011 0.002 0.007 0.003 -0.001 Loadings, Communalities and PSI 0.664 -0.321 0.074 0.550 0.450 0.689 -0.247 -0.193 0.573 0.427 0.493 -0.302 -0.222 0.383 0.617 0.837 0.292 -0.035 0.788 0.212 0.705 0.315 -0.153 0.619 0.381 0.819 0.377 0.105 0.823 0.177 0.661 -0.396 -0.078 0.600 0.400 0.458 -0.296 0.491 0.538 0.462 0.766 -0.427 -0.012 0.769 0.231 Return: The function returns NAG error code or 0 if no error. 70: On entry, parameter matrix had an illegal value. 11: On entry, nfac must not be less than 1: nfac = _value_. On entry, nvar must not be less than 2: nvar = _value_. 17: On entry, m = _value_ while nvar = _value_. These parameters must satisfy m = nvar. On entry, tdx = _value_ while m = _value_. These parameters must satisfy tdx = m. On entry, td. = _value_ while nfac = _value_. These parameters must satisfy td. = nfac. 18: On entry, n = _value_ while nvar = _value_. These parameters must satisfy n > nvar. 19: On entry, nfac = _value_ while nvar = _value_. These parameters must satisfy nfac = nvar. 44: Value _value_ given to eps is not valid. Correct range is machine precision = eps < 1.0. 500: On entry, wt[_value_] = _value_. Constraint: When referenced, all elements of wt must be non-negative. 501: The number of variables, nvar in the analysis = _value_, while number of variables included in the analysis via array isx = _value_. Constraint: these two numbers must be the same. 503: With weighted data, the e.ective number of observations given by the sum of weights = _value_, while the number of variables included in the analysis, nvar = _value_. Constraint: e.ective number of observations > nvar + 1. 427: A singular value decomposition has failed to converge. This is a very unlikely error exit. 284: The conditions for a minimum have not all been satis.ed but a lower point could not be found. Note that in this case all the results are computed. See nag opt bounds 2nd deriv (e04lbc) for further details. 285:The maximum number of iterations, _value_, have been performed. 514: On entry, matrix = Nag DataCorr or matrix = Nag DataCovar and the data matrix is not of full column rank, or matrix = Nag MatCorr Covar and the input correlation/variancecovariance matrix is not positive-de.nite. This exit may also be caused by two of the eigenvalues of S* being equal; this is rare (see Lawley and Maxwell (1971)) and may be due to the data/correlation matrix being almost singular. 73: Memory allocation failed. 74: An internal error has occurred in this function. Check the function call and array sizes. If the function call is correct, please consult NAG for assistance. Additional error messages will output if the optimisation fails to converge or if the options are set incorrectly. Details of these can be found in the nag opt bounds 2nd deriv (e04lbc) document. successfully call of the nag_mv_factor function. */ /** int nag_mv_factor( Nag_FacMat mat, int n, // the number of observations . int m, // the number of variables in the data/correlation/variance-covariance matrix. double x[], // the input matrix. int tdx, // the last dimension of the array x. int nvar, // the number of variables in the factor analysis. int isx[], // indicates whether or not the jth variable is to be included in the factor analysis. int nfac, // the number of factors. double *wtptr, // the weights to be used in the factor analysis. double e[], // the eigenvalues. double stat[], // the test statistics. double com[], // the communalities. double psi[], // the estimates of psi(i), for i = 1, 2, . . . , p. double res[], // the residual correlations. double fl[], // the factor loadings. int tdfl, // the last dimension of the array fl. Nag_E04_Opt *options, double eps, // A lower bound for the value of Øi. ); // Loadings, Communalities and PSI, Eigenvalues, Standardized Communalities, .... */ /** g03ccc computes factor score coefficients from the result of fitting a factor analysis model by maximum likelihood as performed by nag mv factor (g03cac). Example: The example is taken from Lawley and Maxwell (1971). The correlation matrix for 220 observations on six school subjects is input and a factor analysis model with two factors .tted using nag mv factor (g03cac). The factor score coefficients are computed using the regression method. // Note: This one needs the support of the the results of the g03cac Return: The function returns NAG error code or 0 if no error. 70: On entry, parameter method had an illegal value. On entry, parameter rotate had an illegal value. 11: On entry, nfac must not be less than 1: nfac = _value_. 17: On entry, nvar = _value_ while nfac = _value_. These parameters must satisfy nvar = nfac. On entry, td. = _value_ while nfac = _value_. These parameters must satisfy td. = nfac. On entry, tdfs = _value_ while nfac = _value_. These parameters must satisfy tdfs = nfac. 109: On entry, tdr = _value_ while nfac = _value_ and rotate = Nag FacRotate. These parameters must satisfy tdr = nfac when rotate = Nag FacRotate. 55: On entry, e[_value_] = _value_. Constraint: e[_value_] > 1.0. On entry, psi[_value_] = _value_. Constraint: psi[_value_] > 0.0. 73: Memory allocation failed. 74: An internal error has occurred in this function. Check the function call and array sizes. If the function call is correct, please consult NAG for assistance. */ /** int nag_mv_fac_score( Nag_FacScoreMethod method, Nag_FacRotation rotate, int nvar, // the number of observed variables in the factor analysis. int nfac, // the number of factors in the factor analysis. double fl[], // the matrix of unrotated factor loadings, Ë, as returned by nag_mv_factor (g03cac). int tdfl, // the last dimension of the array fl. double psi[], // the diagonal elements of Ø, as returned by nag mv factor (g03cac). double e[], // the eigenvalues of the matrix S*, as returned by nag mv factor (g03cac). double r[], // the orthogonal rotation matrix, R, as returned by nag mv orthomax (g03bac) if rotate = Nag FacRotate. int tdr, // the last dimension of the array r. double fs[], // the matrix of factor score coefficients int tdfs // the last dimension of the array fs. ); // Loadings, Communalities and PSI, Factor score coefficients, Standardized Communalities ... */ /** g03dac computes a test statistic for the equality of within-group covariance matrices and also computes matrices for use in discriminant analysis. Example: The data, taken from Aitchison and Dunsmore (1975), is concerned with the diagnosis of three 'types' of Cushing's syndrome. The variables are the logarithms of the urinary excretion rates (mg/24hr) of two steroid metabolites. Observations for a total of 21 patients are input and the statistics computed by nag mv discrim (g03dac). The printed results show that there is evidence that the within-group variance-covariance matrices are not equal. void test_nag_mv_discrim() { double stat; double det[3], gc[12]; matrix gmean; gmean.SetSize(3,2); matrix x = {{1.1314, 2.4596}, {1.0986, 0.2624}, {0.6419, -2.3026}, {1.3350, -3.2189}, {1.4110, 0.0953}, {0.6419, -0.9163}, {2.1163, 0.0000}, {1.3350, -1.6094}, {1.3610, -0.5108}, {2.0541, 0.1823}, {2.2083, -0.5108}, {2.7344, 1.2809}, {2.0412, 0.4700}, {1.8718, -0.9163}, {1.7405, -0.9163}, {2.6101, 0.4700}, {2.3224, 1.8563}, {2.2192, 2.0669}, {2.2618, 1.1314}, {3.9853, 0.9163}, {2.7600, 2.0281}}; double wt[21]; double df; double *wtptr = NULL; double sig; int nvar = 2; int i, j, m = 2, n = 21; int isx[2], nig[3]; int ing[21] = {1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3}; int ng = 3; int tdgmean = 2; int tdx = 2; char weight = 'U'; for (j = 0; j < m; ++j) isx[j] = 1; nag_mv_discrim(n, m, x, tdx, isx, nvar, ing, ng, wtptr, nig, gmean, tdgmean, det, gc, &stat, &df, &sig); printf("\nGroup means\n\n"); for (i = 0; i < ng; ++i) { for (j = 0; j < nvar; ++j) printf("%10.4f",gmean[i][j]); printf("\n"); } printf("\nLOG of determinants\n\n"); for (j = 0; j < ng; ++j) printf("%10.4f",det[j]); printf("\n\n%s%7.4f\n","stat = ",stat); printf("%s%7.4f\n"," df = ",df); printf("%s%7.4f\n"," sig = ",sig); } The output is following: Group means 1.0433 -0.6034 2.0073 -0.2060 2.7097 1.5998 LOG of determinants -0.8273 -3.0460 -2.2877 stat = 19.2410 df = 6.0000 sig = 0.0038 Return: The function returns NAG error code or 0 if no error. 11: On entry, n must not be less than 1: n = _value_. On entry, nvar must not be less than 1: nvar = _value_. On entry, ng must not be less than 2: ng = _value_. 17: On entry, m = _value_ while nvar = _value_. These parameters must satisfy m = nvar. On entry, tdx = _value_ while m = _value_. These parameters must satisfy tdx = m. On entry, tdg = _value_ while nvar = _value_. These parameters must satisfy tdg = nvar. 112: On entry, ing[_value_] = _value_, ng = _value_. Constraint: 1 = ing[i - 1] = ng, i = 1, 2, . . ., n. 500: On entry, wt[_value_] = _value_. Constraint: when referenced, all elements of wt must be non-negative. 501: The number of variables, nvar in the analysis = _value_, while number of variables included in the analysis via array isx = _value_. Constraint: these two numbers must be the same. 515: On entry, group _value_ has _value_ e.ective observations. Constraint: in each group the e.ective number of observations must be = 1. 516: On entry, group _value_ has _value_ members, while nvar = _value_. Constraint: number of members in each group = nvar 520: The variables in group _value_ are not of full rank. 519: The variables are not of full rank. 73: Memory allocation failed. 74: An internal error has occurred in this function. Check the function call and array sizes. If the function call is correct, please consult NAG for assistance. successfully call of the nag_mv_discrim function. */ int nag_mv_discrim( int n, // the number of observations. int m, // the number of variables in the data array x. double x[], // contain the kth observation for the lth variable, for k = 1, 2, . . ., n; l = 1, 2, . . ., m . int tdx, // the last dimension of the array x. int isx[], // indicates whether or not the lth variable in x is to be included in the variance-covariance matrices. int nvar, // the number of variables in the variance-covariance matrices. int ing[], // indicates to which group the kth observation belongs, for k = 1, 2, . . ., n. int ng, // the number of groups. double *wtptr, // the weights to be used in the analysis and the effective number of observations for a group is the sum of the weights of the observations in that group. int nig[], // the number of observations in the jth group, for j = 1, 2, . . . , ng double gmean[], // the jth row of gmean contains the means of the p selected variables for the jth group, for j = 1, 2, . . . , ng. int tdg, // the last dimension of the array gmean. double det[], // the logarithm of the determinants of the within-group variance-covariance matrices. double gc[], // the first p(p+1)/2 elements of gc contain R and the remaining ng blocks of p(p+1)/2 elements contain the Rj matrices. double *stat, // the likelihood-ratio test statistic, G. double *df, // the degrees of freedom for the distribution of G. double *sig // the significance level for G. ); /** g03dbc computes Mahalanobis squared distances for groups or pooled variance-covariance matrices. It is intended for use after nag mv discrim (g03dac). Example: The data, taken from Aitchison and Dunsmore (1975), are concerned with the diagnosis of three 'types' of Cushing's syndrome. The variables are the logarithms of the urinary excretion rates (mg/24hr) of two steroid metabolites. Observations for a total of 21 patients are input and the group means and R matrices are computed by nag mv discrim (g03dac). A further six observations of unknown type are input, and the distances from the group means of the 21 patients of known type are computed under the assumption that the within-group variance-covariance matrices are not equal. These results are printed and indicate that the .rst four are close to one of the groups while observations 5 and 6 are some distance from any group. void test_nag_mv_discrim_mahaldist() { matrix d, gmean; d.SetSize(21,3); gmean.SetSize(3,2); matrix x = {{1.1314, 2.4596}, {1.0986, 0.2624}, {0.6419, -2.3026}, {1.3350, -3.2189}, {1.4110, 0.0953}, {0.6419, -0.9163}, {2.1163, 0.0000}, {1.3350, -1.6094}, {1.3610, -0.5108}, {2.0541, 0.1823}, {2.2083, -0.5108}, {2.7344, 1.2809}, {2.0412, 0.4700}, {1.8718, -0.9163}, {1.7405, -0.9163}, {2.6101, 0.4700}, {2.3224, 1.8563}, {2.2192, 2.0669}, {2.2618, 1.1314}, {3.9853, 0.9163}, {2.7600, 2.0281}}; double det[3], gc[12]; double wt[21]; double stat; double df; double sig; double *wtptr; int nobs, nvar; int isx[2], nig[3]; int ing[21] = {1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3}; int i, j, m = 2, n = 21; int ng = 3; int tdd = 3; int tdgmean = 2; int tdx = 2; char char_equal; char weight = 'U'; Nag_GroupCovars equal; wtptr = NULL; nvar = 2; for (j = 0; j < m; ++j) isx[j] = 1; nag_mv_discrim(n, m, x, tdx, isx, nvar, ing, ng, wtptr, nig, gmean, tdgmean, det, gc, &stat, &df, &sig); nobs = 6; char_equal = 'U'; x[0][0] = 1.6292; x[0][1] = -0.9163; x[1][0] = 2.5572; x[1][1] = 1.6094; x[2][0] = 2.5649; x[2][1] = -0.2231; x[3][0] = 0.9555; x[3][1] = -2.3026; x[4][0] = 3.4012; x[4][1] = -2.3026; x[5][0] = 3.0204; x[5][1] = -0.2231; if(char_equal == 'E') equal = Nag_EqualCovar; else if (char_equal == 'U') equal = Nag_NotEqualCovar; nag_mv_discrim_mahaldist(equal, Nag_SamplePoints, nvar, ng, gmean, tdgmean, gc, nobs, m, isx, x, tdx, d, tdd); printf("\n Obs Distances\n\n"); for (i = 0; i < nobs; ++i) { printf(" %3ld",i+1); for (j = 0; j < ng; ++j) printf("%10.3f",d[i][j]); printf("\n"); } } The Output is following: Obs Distances 1 3.339 0.752 50.928 2 20.777 5.656 0.060 3 21.363 4.841 19.498 4 0.718 6.280 124.732 5 55.000 88.860 71.785 6 36.170 15.785 15.749 Return: The function returns NAG error code or 0 if no error. 70: On entry, parameter equal had an illegal value. On entry, parameter mode had an illegal value. 11: On entry, nvar must not be less than 1: nvar = _value_. On entry, ng must not be less than 2: ng = _value_. 17: On entry, tdg = _value_ while nvar = _value_. These parameters must satisfy tdg = nvar. On entry, tdd = _value_ while ng = _value_. These parameters must satisfy tdd = ng. 108: On entry, nobs = _value_ while mode = Nag SamplePoints. These parameters must satisfy nobs = 1 w hen mode = Nag SamplePoints. 109: On entry, m = _value_ while nvar = _value_ and mode = Nag SamplePoints. These parameters must satisfy m = nvar when mode = Nag SamplePoints. On entry, tdx = _value_ while m = _value_ and mode = Nag SamplePoints. These parameters must satisfy tdx = max(1,m) w hen mode = Nag SamplePoints. 502: The number of variables, nvar in the analysis = _value_, while number of variables included in the analysis via array isx = _value_. Constraint: These two numbers must be the same when mode = Nag SamplePoints. 517: A diagonal element of R is zero when equal = Nag EqualCovar. 518: A diagonal element of R is zero for some j, w hen equal = Nag NotEqualCovar. 73: Memory allocation failed. 74: An internal error has occurred in this function. Check the function call and array sizes. If the function call is correct, please consult NAG for assistance. successfully call of the nag_mv_discrim_mahaldist function. */ int nag_mv_discrim_mahaldist( Nag_GroupCovars equal, Nag_MahalDist mode, int nvar, // the number of variables, p, in the variance-covariance matrices as specified to nag mv discrim (g03dac). int ng, // the number of groups. double gmean[], // the jth row of gmean contains the means of the p selected variables for the jth group, for j = 1, 2, . . . , ng. int tdg, // the last dimension of the array gmean. double gc[], // the first p(p + 1)/2 elements of gc should contain the upper triangular matrix R and the next ng blocks of p(p+1)/2 elements should contain the upper triangular matrices Rj. int nobs, // if mode = Nag SamplePoints the number of sample points in x for which distances are to be calculated. int m, // if mode = Nag SamplePoints the number of variables in the data array x. int isx[], // if mode = Nag SamplePoints, isx[l-1] indicates if the lth variable in x is to be included in the distance calculations. double x[], // contain the kth sample value for the lth variable for k = 1, 2, . . . ,nobs; l = 1, 2, . . . ,m. int tdx, // the last dimension of the array x. double d[], // the squared distances. int tdd // the last dimension of the array dd. ); /** g03dcc allocates observations to groups according to selected rules. It is intended for use after nag mv discrim (g03dac). Example The data, taken from Aitchison and Dunsmore (1975), are concerned with the diagnosis of three 'types' of Cushing's syndrome. The variables are the logarithms of the urinary excretion rates (mg/24hr) of two steroid metabolites. Observations for a total of 21 patients are input and the group means and R matrices are computed by nag mv discrim (g03dac). A further six observations of unknown type are input and allocations made using the predictive approach and under the assumption that the within-group covariance matrices are not equal. The posterior probabilities of group membership, q[j], and the atypicality index are printed along with the allocated group. The atypicality index shows that observations 5 and 6 do not seem to be typical of the three types present in the initial 21 observations. void test_nag_mv_discrim_group() { double stat; matrix ati, gmean, p; ati.SetSize(21,3); gmean.SetSize(3,2); p.SetSize(21,3); matrix x = {{1.1314, 2.4596}, {1.0986, 0.2624}, {0.6419, -2.3026}, {1.3350, -3.2189}, {1.4110, 0.0953}, {0.6419, -0.9163}, {2.1163, 0.0000}, {1.3350, -1.6094}, {1.3610, -0.5108}, {2.0541, 0.1823}, {2.2083, -0.5108}, {2.7344, 1.2809}, {2.0412, 0.4700}, {1.8718, -0.9163}, {1.7405, -0.9163}, {2.6101, 0.4700}, {2.3224, 1.8563}, {2.2192, 2.0669}, {2.2618, 1.1314}, {3.9853, 0.9163}, {2.7600, 2.0281}}; double det[3], gc[12]; double prior[3], wt[21]; double df; double sig; double *wtptr = NULL; int nobs, nvar; int i, j, m, n; int iag[21], isx[2], nig[3]; int ing[21] = {1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3}; int ng = 3; int tdgmean = 2; int tdp = 3; int tdx = 2; Boolean atiq = true; char char_type; char char_equal; char weight = 'U'; Nag_DiscrimMethod type; Nag_GroupCovars equal; n = 21; m = 2; nvar = 2; for (j = 0; j < m; ++j) isx[j] = 1; nag_mv_discrim(n, m, x, tdx, isx, nvar, ing, ng, wtptr, nig, gmean, tdgmean, det, gc, &stat, &df, &sig); nobs = 6; char_equal ='U'; char_type = 'P'; x[0][0] = 1.6292; x[0][1] = -0.9163; x[1][0] = 2.5572; x[1][1] = 1.6094; x[2][0] = 2.5649; x[2][1] = -0.2231; x[3][0] = 0.9555; x[3][1] = -2.3026; x[4][0] = 3.4012; x[4][1] = -2.3026; x[5][0] = 3.0204; x[5][1] = -0.2231; if (char_type == 'E') type = Nag_DiscrimEstimate; if(char_type == 'P') type = Nag_DiscrimPredict; if(char_equal == 'E') equal = Nag_EqualCovar; if (char_equal == 'U') equal = Nag_NotEqualCovar; nag_mv_discrim_group(type, equal, Nag_EqualPrior, nvar, ng, nig, gmean, tdgmean, gc, det, nobs, m, isx, x, tdx, prior, p, tdp, iag, atiq, ati); printf("\n"); printf(" Obs Posterior Allocated "); printf(" Atypicality "); printf("\n"); printf(" probabilities to group index "); printf("\n"); printf("\n"); for (i = 0; i < nobs; ++i) { printf(" %6ld ",i+1); for (j = 0; j < ng; ++j) printf("%6.3f",p[i][j]); printf(" %6ld ",iag[i]); for (j = 0; j < ng; ++j) printf("%6.3f",ati[i][j]); printf("\n"); } } The output is following: Obs Posterior Allocated Atypicality probabilities to group index 1 0.094 0.905 0.002 2 0.596 0.254 0.975 2 0.005 0.168 0.827 3 0.952 0.836 0.018 3 0.019 0.920 0.062 2 0.954 0.797 0.912 4 0.697 0.303 0.000 1 0.207 0.860 0.993 5 0.317 0.013 0.670 3 0.991 1.000 0.984 6 0.032 0.366 0.601 3 0.981 0.978 0.887 Return: The function returns NAG error code or 0 if no error. 70: On entry, parameter type had an illegal value. On entry, parameter equal had an illegal value. On entry, parameter priors had an illegal value. 11: On entry, nvar must not be less than 1: nvar = _value_. On entry, ng must not be less than 2: ng = _value_. On entry, nobs must not be less than 1: nobs = _value_. 17: On entry, m = _value_ while nvar = _value_. These parameters must satisfy m = nvar. On entry, tdx = _value_ while m = _value_. These parameters must satisfy tdx = m. On entry, tdp = _value_ while ng = _value_. These parameters must satisfy tdp = ng. On entry, tdg = _value_ while nvar = _value_. These parameters must satisfy tdg = nvar. 501: The number of variables, nvar in the analysis = _value_, while number of variables included in the analysis via array isx = _value_. Constraint: these two numbers must be the same. 111: On entry, nig[_value_] = _value_. Constraint: nig[i - 1] > 0, i = 1, 2, . . . ,ng when equal = Nag EqualCovar. 112: On entry, nig[_value_] = _value_, nvar = _value_. Constraint: nig[i - 1] > nvar, i = 1, 2, . . . ,ng when equal = Nag NotEqualCovar. 115: On entry, prior[_value_] = _value_. Constraint: prior[j - 1] > 0, j=1,2,. . . ,ng when priors = Nag UserPrior. 522: On entry,summation of prior[j - 1] = _value_. Constraint: summation of prior[j - 1] must be within 10 × machine precision of 1 when priors = Nag UserPrior. 521: On entry, summation of nig[j - 1] = _value_, ng = _value_, nvar = _value_. Constraint: summation of nig[j - 1] > ng + nvar when equal = Nag EqualCovar. 517: A diagonal element of R is zero when equal = Nag EqualCovar. 518: A diagonal element of R is zero for some j, when equal = Nag NotEqualCovar 73: Memory allocation failed. 74: An internal error has occurred in this function. Check the function call and array sizes. If the function call is correct, please consult NAG for assistance. successfully call of the nag_mv_discrim_group function. */ int nag_mv_discrim_group( Nag_DiscrimMethod type, Nag_GroupCovars equal, Nag_PriorProbability priors, int nvar, // the number of variables, p, in the variance-covariance matrices as specified to nag mv discrim (g03dac). int ng, // the number of groups. int nig[], // the number of observations in each group training set. double gmean[], // the jth row of gmean contains the means of the p variables for the jth group, for j = 1, 2,.., nj. int tdg, // the last dimension of the array gmean. double gc[], // the first p(p + 1)/2 elements of gc should contain the upper triangular matrix R and the next ng blocks of p(p + 1)/2 elements should contain the upper triangular matrices Rj . double det[], // if equal = Nag NotEqualCovar the logarithms of the determinants of the within-group variance-covariance matrices as returned by nag mv discrim (g03dac). int nobs, // the number of observations in x which are to be allocated. int m, // the number of variables in the data array x. int isx[], // indicates if the lth variable in x is to be included in the distance calculations. double x[], // contain the kth observation for the lth variable, for k = 1, 2, . . . ,nobs; l = 1, 2, . . . ,m. int tdx, // the last dimension of the array x. double prior[], // if priors = Nag UserPrior, the prior probabilities for the ng groups. double p[], // contains the posterior probability p[k][j] for allocating the kth observation to the jth group, for k = 1, 2, . . . ,nobs; j = 1, 2, . . ., ng. int tdp, // the last dimension of the array p and ati. int iag[], // the groups to which the observations have been allocated. Boolean atiq, // atiq must be TRUE if atypicality indices are required. . double ati[] // the atypicality index for the kth observation with respect to the jth group, // for k = 1, 2, . . . ,nobs; j = 1, 2, . . ., ng if atiq is TRUE, ); /** g03eac computes a distance (dissimilarity) matrix. Example: A data matrix of five observations and three variables is read in and a distance matrix is calculated from variables 2 and 3 using squared Euclidean distance with no scaling. This matrix is then printed. void test_nag_mv_distance_mat() { double d[45], s[10]; matrix x; x.SetSize(10,10); int isx[10]; int i, j, m, n; int tdx = 10; char char_scale; char char_update; char char_dist; Nag_MatUpdate update; Nag_DistanceType dist; Nag_VarScaleType scale; n = 5; m = 3; char_update = 'I'; char_dist = 'S'; char_scale = 'U'; x[0][0] = 1.0; x[0][1] = 1.0; x[0][2] = 1.0; x[1][0] = 2.0; x[1][1] = 1.0; x[1][2] = 2.0; x[2][0] = 3.0; x[2][1] = 6.0; x[2][2] = 3.0; x[3][0] = 4.0; x[3][1] = 8.0; x[3][2] = 2.0; x[4][0] = 5.0; x[4][1] = 8.0; x[4][2] = 0.0; for (i = 0; i < m; ++i) isx[i] = 1; isx[0] = 0; for (i = 0; i < m; ++i) s[i] = 1.0; if (char_update == 'U') update = Nag_MatUp; if (char_update == 'I') update = Nag_NoMatUp; if (char_dist == 'A') dist = Nag_DistAbs; if (char_dist == 'E') dist = Nag_DistEuclid; if (char_dist == 'S') dist = Nag_DistSquared; if (char_scale == 'S') scale = Nag_VarScaleStd; if (char_scale == 'R') scale = Nag_VarScaleRange; if (char_scale == 'G') scale = Nag_VarScaleUser; if (char_scale == 'U') scale = Nag_NoVarScale; nag_mv_distance_mat(update, dist, scale, n, m, x, tdx, isx, s, d); printf("\n"); printf(" Distance matrix "); printf("\n"); printf("\n"); printf(" %s\n"," 1 2 3 4"); printf("\n"); for (i = 2; i <= n; ++i) { printf("%2ld ",i); for (j=(i-1)*(i-2)/2+1; j<=i*(i - 1)/2; ++j) printf("%5.2f ",d[j-1]); printf("\n"); } } The output is following: Distance matrix 1 2 3 4 2 1.00 3 29.00 26.00 4 50.00 49.00 5.00 5 50.00 53.00 13.00 4.00 Return: The function returns NAG error code or 0 if no error. 70: On entry, parameter dist had an illegal value. On entry, parameter update had an illegal value. On entry, parameter scale had an illegal value. 11: On entry, n must not be less than 2: n = _value_. 12: On entry, m must not be less than or equal to 0: m = _value_. 17: On entry, tdx = _value_ while m = _value_. These parameters must satisfy tdx = m. 111: On entry, isx[_value_] = _value_. Constraint: isx[i - 1] > 0, for at least one i, i = 1, 2, . . . ,m. 115: On entry, d[_value_] = _value_. Constraint: d[i - 1] = 0, i = 1, 2, . . . ,n*(n-1)/2 when update = Nag MatUp. On entry, s[_value_] = _value_. Constraint: s[j - 1] > 0, j = 1, 2, . . . ,m when scale = Nag VarScaleUser and isx[j - 1] > 0. 523: On entry, scale = Nag VarScaleRange or scale = Nag VarScaleStd, and x[i - 1][j - 1] = x[i][j - 1], for i = 1, 2, . . . , n - 1, for some j with isx[i - 1] > 0. 74: An internal error has occurred in this function. Check the function call and array sizes. If the function call is correct, please consult NAG for assistance. successfully call of the nag_mv_distance_mat function. */ int nag_mv_distance_mat( Nag_MatUpdate update, Nag_DistanceType dist, Nag_VarScaleType scale, int n, // the number of observations. int m, // the total number of variables in array x. double x[], // the value of the jth variable for the ith object, for i = 1, 2, . . . , n; j = 1, 2, . . . ,m. int tdx, // the last dimension of the array x. int isx[], // indicates whether or not the jth variable in x is to be included in the distance computations. double s[], // Input:if scale = Nag VarScaleUser and isx[j - 1] > 0 then s[j - 1] must contain the scaling for variable j, for j = 1, 2, . . . ,m. Output:if scale = Nag VarScaleStd and isx[j - 1] > 0 then s[j - 1] contains the standard deviation of the variable in the jth column of x. double d[] // Input:if update = Nag MatUp then d must contain the strictly lower triangle of the distance matrix D to be updated. Output:the strictly lower triangle of the distance matrix D. ); /** g03ecc performs hierarchical cluster analysis. Example: Data consisting of three variables on five objects are read in. Euclidean squared distances based on two variables are computed using nag mv distance mat (g03eac), the objects are clustered using nag mv hierar cluster analysis and the dendrogram computed using nag mv dendrogram (g03ehc). The dendrogram is then printed. void test_nag_mv_hierar_cluster_analysis() { double cd[9], d[45], dord[10], s[10]; matrix x; x.SetSize(10,10); double dmin_, dstep, ydist; int ilc[9], iord[10], isx[10], iuc[9]; int nsym; int i, j; int m = 3, n = 5; int int_method; int tdx = 10; char name[10] = "ABCDE"; char char_dist; char **c= NULL; char char_scale; char char_update; Nag_ClusterMethod method; Nag_MatUpdate update; Nag_DistanceType dist; Nag_VarScaleType scale; int_method = 5; if (int_method == 1) method = Nag_SingleLink; else if (int_method == 2) method = Nag_CompleteLink; else if (int_method == 3) method = Nag_GroupAverage; else if (int_method == 4) method = Nag_Centroid; else if (int_method == 5) method = Nag_Median; else method = Nag_MinVariance; char_update = 'I'; if (char_update == 'U') update = Nag_MatUp; else update = Nag_NoMatUp; char_dist = 'S'; if (char_dist == 'A') dist = Nag_DistAbs; else if (char_dist == 'E') dist = Nag_DistEuclid; else dist = Nag_DistSquared; char_scale = 'U'; if (char_scale == 'S') scale = Nag_VarScaleStd; else if (char_scale == 'R') scale = Nag_VarScaleRange; else if(char_scale == 'G') scale = Nag_VarScaleUser; else scale = Nag_NoVarScale; x[0][0] = 1; x[0][1] = 5.0; x[0][2] = 2.0; x[1][0] = 2; x[1][1] = 1.0; x[1][2] = 1.0; x[2][0] = 3; x[2][1] = 4.0; x[2][2] = 3.0; x[3][0] = 4; x[3][1] = 1.0; x[3][2] = 2.0; x[4][0] = 5; x[4][1] = 5.0; x[4][2] = 0.0; for (i = 0; i < m; ++i) isx[i] = 1; isx[0] = 0; for (i = 0; i < m; ++i) s[i] = 1.0; nag_mv_distance_mat(update, dist, scale, n, m, x, tdx, isx, s, d); nag_mv_hierar_cluster_analysis(method, n, d, ilc, iuc, cd, iord, dord); printf("\n Distance Clusters Joined\n\n"); for (i = 1; i <= n-1; ++i) printf("%10.3f %3c%3c\n",cd[i-1], name[ilc[i-1]-1], name[iuc[i-1]-1]); nsym = 20; dmin_ = 0.0; dstep = cd[n - 2] / (1.0 * nsym); nag_mv_dendrogram(Nag_DendSouth, n, dord, dmin_, dstep, nsym, &c); printf("\n"); printf("Dendrogram "); printf("\n"); printf("\n"); ydist = cd[n - 2]; for (i = 0; i < nsym; ++i) { if ((i+1) % 3 == 1) { printf("%10.3f%6s",ydist,""); printf("%s",c[i]); printf("\n"); } else { printf("%16s%s","", c[i]); printf("\n"); } ydist -= dstep; } printf("\n"); printf("%14s",""); for (i = 0; i < n; ++i) printf("%3c",name[iord[i]-1]); printf("\n"); nag_mv_dend_free(&c); } The output is as following: Distance Clusters Joined 1.000 B D 2.000 A C 6.500 A E 14.125 A B Dendrogram 14.125 ------- I I I I 12.006 I I I I I I 9.887 I I I I I I 7.769 I I ---* I I I I 5.650 I I I I I I I I I 3.531 I I I I I I ---* I I 1.412 I I I ---* I I I I I A C E B D Return: The function returns NAG error code or 0 if no error. 70: On entry, parameter method had an illegal value. 11: On entry, n must not be less than 2: n = _value_. 115: On entry, d[_value_] = _value_. Constraint: d[i - 1] = 0.0, i = 1, 2, . . . , n * (n - 1)/2. 524: A true dendrogram cannot be formed because the distances at which clusters have merged are not increasing for all steps, i.e., cd[i - 1] < cd[i - 2] for some i = 2, 3, . . . , n - 1. This can occur for the Nag Centroid and Nag Median methods. 73: Memory allocation failed. 74: An internal error has occurred in this function. Check the function call and array sizes. If the function call is correct, please consult NAG for assistance. successfully call of the nag_mv_hierar_cluster_analysis function. */ int nag_mv_hierar_cluster_analysis( Nag_ClusterMethod method, int n, // the number of objects. double d[], // the strictly lower triangle of the distance matrix. int ilc[], // contains the number, j, of the cluster merged with cluster k (see iuc), j < k, at step l for l = 1, 2, . n - 1. int iuc[], // contains the number, k, of the cluster merged with cluster j, j < k, at step l for l = 1, 2, . . . , n - 1. double cd[], // the distance djk, between clusters j and k, j < k, merged at step l for l = 1, 2, . . . , n - 1. int iord[], // the objects in dendrogram order. double dord[] // the clustering distances corresponding to the order in iord. ); /** g03efc performs K-means cluster analysis. Example: The data consists of observations of .ve variables on twenty soils (Kendall and Stuart (1976)). The data is read in, the K-means clustering performed and the results printed. void test_nag_mv_kmeans_cluster_analysis() { double css[3], csw[3], wt[20]; double *wtptr; int nvar, i, j, k; int m, n; int inc[20], isx[5], nic[5]; int maxit; int tdc = 5; int tdx = 5; char weight; matrix x = {{77.3, 13.0, 9.7, 1.5, 6.4}, {82.5, 10.0, 7.5, 1.5, 6.5}, {66.9, 20.6, 12.5, 2.3, 7.0}, {47.2, 33.8, 19.0, 2.8, 5.8}, {65.3, 20.5, 14.2, 1.9, 6.9}, {83.3, 10.0, 6.7, 2.2, 7.0}, {81.6, 12.7, 5.7, 2.9, 6.7}, {47.8, 36.5, 15.7, 2.3, 7.2}, {48.6, 37.1, 14.3, 2.1, 7.2}, {61.6, 25.5, 12.9, 1.9, 7.3}, {58.6, 26.5, 14.9, 2.4, 6.7}, {69.3, 22.3, 8.4, 4.0, 7.0}, {61.8, 30.8, 7.4, 2.7, 6.4}, {67.7, 25.3, 7.0, 4.8, 7.3}, {57.2, 31.2, 11.6, 2.4, 6.5}, {67.2, 22.7, 10.1, 3.3, 6.2}, {59.2, 31.2, 9.6, 2.4, 6.0}, {80.2, 13.2, 6.6, 2.0, 5.8}, {82.2, 11.1, 6.7, 2.2, 7.2}, {69.7, 20.7, 9.6, 3.1, 5.9}}; matrix cmeans = {{82.5, 10.0, 7.5, 1.5, 6.5}, {47.8, 36.5, 15.7, 2.3, 7.2}, {67.2, 22.7, 10.1, 3.3, 6.2}}; weight = 'U'; n = 20; m = 5; nvar = 5; k = 3; maxit = 10; for (j = 0; j < m; ++j) isx[j] = 1; nag_mv_kmeans_cluster_analysis(n, m, x, tdx, isx, nvar, k, cmeans, tdc, wtptr, inc, nic, css, csw, maxit); printf("\nThe cluster each point belongs to\n"); for (i = 0; i < n; ++i) { printf(" %6ld",inc[i]); if((i + 1) % 10 == 0) printf("\n"); } printf("\n\nThe number of points in each cluster\n"); for (i = 0; i < k; ++i) printf(" %6ld",nic[i]); printf("\n\nThe within-cluster sum of weights of each cluster\n"); for (i = 0; i < k; ++i) printf(" %9.2f",csw[i]); printf("\n\nThe within-cluster sum of squares of each cluster\n\n"); for (i = 0; i < k; ++i) printf(" %13.4f",css[i]); printf("\n\nThe final cluster centres\n"); printf(" 1 2 3 4 5\n"); for (i = 0; i < k; ++i) { printf(" %5ld ",i+1); for (j = 0; j < nvar; ++j) printf("%8.4f",cmeans[i][j]); printf("\n"); } } The output is following: The cluster each point belongs to 1 1 3 2 3 1 1 2 2 3 3 3 3 3 3 3 3 1 1 3 The number of points in each cluster 6 3 11 The within-cluster sum of weights of each cluster 6.00 3.00 11.00 The within-cluster sum of squares of each cluster 46.5717 20.3800 468.8964 The final cluster centres 1 2 3 4 5 1 81.1833 11.6667 7.1500 2.0500 6.6000 2 47.8667 35.8000 16.3333 2.4000 6.7333 3 64.0455 25.2091 10.7455 2.8364 6.6545 Return: The function returns NAG error code or 0 if no error. 11: On entry, n must not be less than 2: n = _value_. On entry, k must not be less than 2: k = _value_. On entry, nvar must not be less than 1: nvar = _value_. 12: On entry, maxit must not be less than or equal to 0: maxit = _value_. 17: On entry, m = _value_ while nvar = _value_. These parameters must satisfy m = nvar. On entry, tdx = _value_ while m = _value_. These parameters must satisfy tdx = m. On entry, tdc = _value_ while nvar = _value_. These parameters must satisfy tdc = nvar. 501: The number of variables, nvar in the analysis = _value_, while number of variables included in the analysis via array isx = _value_. Constraint: these two numbers must be the same. 500: On entry, wt[_value_] = _value_. Constraint: When referenced, all elements of wt must be non-negative. 525: At least two elements of wt must be greater than zero. 526: At least one cluster is empty after the initial assignment. Try a di.erent set of initial cluster centres in cmeans and also consider decreasing the value of k. The empty clusters may be found by examining the values in nic. 3: Too many iterations (_value_). Convergence has not been achieved within the maximum number of iterations given by maxit. Try increasing maxit and, if possible, use the returned values in cmeans as the initial cluster centres. 73: Memory allocation failed. 74: An internal error has occurred in this function. Check the function call and array sizes. If the function call is correct, please consult NAG for assistance. successfully call of the nag_mv_kmeans_cluster_analysis function. */ int nag_mv_kmeans_cluster_analysis( int n, // the number of observations. int m, // the number of variables in the array x. double x[], // the value of jth variable for the ith object for i = 1, 2, . . ., n; j = 1, 2, . . . ,m. int tdx, // the last dimension of the array x. int isx[], // indicates whether or not the jth variable is to be included in the analysis. int nvar, // the number of variables included in the sum of squares calculations. int k, // the number of clusters. double cmeans[], // Input: contain the value of the jth variable for the ith initial cluster centre, for i = 1, 2, . . .,K; j = 1, 2, . . . , p. Output: contains the value of the jth variable for the ith computed cluster centre, for i = 1, 2, . . .,K; j = 1, 2, . . . , p. int tdc, // the last dimension of the array cmeans. double *wtptr, // contain the weights to be used in the analysis. int inc[], // the cluster to which the ith object has been allocated, for i = 1, 2, . . . , n. int nic[], // the number of objects in the ith cluster, for i = 1, 2, . . .,K. css[k] double css[], // the within-cluster (weighted) sum of squares of the ith cluster, for i = 1, 2, . . .,K. double csw[], // s the within-cluster sum of weights of the ith cluster, for i = 1, 2, . . .,K. int maxit // the maximum number of iterations allowed in the analysis. ); /** g03ehc produces a dendrogramfrom the results of nag mv hierar cluster analysis. Example: Data consisting of three variables on five objects are read in. Euclidean squared distances are computed using nag mv distance mat (g03eac) and median clustering performed by nag mv hierar cluster analysis (g03ecc). nag mv dendrogram (g03ehc) is used to produce a dendrogram with orientation east and a dendrogram with orientation south. The two dendrograms are printed. Note the use of nag mv dend free (g03xzc)to free the memory allocated internally to the character array pointed to by c. void test_nag_mv_dendrogram() { double cd[9], d[45], dord[10], s[10]; matrix x; x.SetSize(10,10); double dmin_; double dstep; int ilc[9], iord[10], isx[10], iuc[9]; int nsym; int i, j; int m = 3, n = 5; int int_method; int tdx; char char_dist; char char_scale; char char_update; char **c = NULL; Nag_ClusterMethod method; Nag_MatUpdate update; Nag_DistanceType dist; Nag_VarScaleType scale; tdx = 10; x[0][0] = 1.0; x[0][1] = 1.0; x[0][2] = 1.0; x[1][0] = 2.0; x[1][1] = 1.0; x[1][2] = 2.0; x[2][0] = 3.0; x[2][1] = 6.0; x[2][2] = 3.0; x[3][0] = 4.0; x[3][1] = 8.0; x[3][2] = 2.0; x[4][0] = 5.0; x[4][1] = 8.0; x[4][2] = 0.0; for (i = 0; i < m; ++i) isx[i] = 1; isx[0] = 0; for (i = 0; i < m; ++i) s[i] = 1; int_method = 5; if (int_method == 1) method = Nag_SingleLink; else if (int_method == 2) method = Nag_CompleteLink; else if (int_method == 3) method = Nag_GroupAverage; else if (int_method == 4) method = Nag_Centroid; else if (int_method == 5) method = Nag_Median; else method = Nag_MinVariance; char_update = 'I'; char_dist = 'S'; char_scale = 'U'; if (char_update == 'U') update = Nag_MatUp; else update = Nag_NoMatUp; if (char_dist == 'A') dist = Nag_DistAbs; else if (char_dist == 'E') dist = Nag_DistEuclid; else dist = Nag_DistSquared; if (char_scale == 'S') scale = Nag_VarScaleStd; else if (char_scale == 'R') scale = Nag_VarScaleRange; else if (char_scale == 'G') scale = Nag_VarScaleUser; else scale = Nag_NoVarScale; dmin_ = 0.0; dstep = 1.1; nsym = 40; nag_mv_distance_mat(update, dist, scale, n, m, x, tdx, isx, s, d); nag_mv_hierar_cluster_analysis(method, n, d, ilc, iuc, cd, iord, dord); nag_mv_dendrogram(Nag_DendEast, n, dord, dmin_, dstep, nsym, &c); printf("\nDendrogram, Orientation East\n\n"); for (i=0; i< n; i++) printf("%s\n", c[i]); dmin_ = 0.0; dstep = 1.0; nsym = 40; nag_mv_dend_free(&c); nag_mv_dendrogram(Nag_DendSouth, n, dord, dmin_, dstep, nsym, &c); printf("\n\n Dendrogram, Orientation South\n\n"); for (i=0; i< nsym; i++) { printf("%s\n", c[i]); } nag_mv_dend_free(&c); } The output is as following: Dendrogram, Orientation East ...............................( ( ....... ( ( ... (........................(...(... Dendrogram, Orientation South ---------- I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I ------* I I I I I I I I I I I ---* I I I I I I I I ---* I I I I I I I I Return: The function returns NAG error code or 0 if no error. 70: On entry, parameter orient had an illegal value. 11: On entry, n must not be less than 2: n = _value_. On entry, nsym must not be less than 1: nsym = _value_. 5: On entry, dmin must not be less than 0.0: dmin = _value_. 6: On entry, dstep must not be less than or equal to 0.0: dstep = _value_. 527: On entry, n = _value_, dord[_value_] = _value_. Constraint: dord[n-1] = dord[i - 1], i = 1, 2, . . . ,n-1. 74: An internal error has occurred in this function. Check the function call and array sizes. If the function call is correct, please consult NAG for assistance. successfully call of the nag_mv_dendrogram function. */ int nag_mv_dendrogram( Nag_DendOrient orient, int n, // the number of objects in the cluster analysis. double dord[], // contains the distances, in dendrogram order, at which clustering takes place. double dmin_, // the clustering distance at which the dendrogram begins. double dstep, // the distance represented by one symbol of the dendrogram. int nsym, // the number of character positions used in the dendrogram. Hence the clustering distance at which the dendrogram terminates is given by dmin + nsym × dstep. char ***c // a pointer to an array of character pointers, containing consecutive lines of the dendrogram. ); /** g03ejc computes a cluster indicator variable from the results of nag mv hierar cluster analysis (g03ecc). Example: Data consisting of three variables on five objects are input. Euclidean squared distances are computed using nag mv distance mat (g03eac) and median clustering performed using nag mv hierar cluster analysis (g03ecc). A dendrogram is produced by nag mv dendrogram (g03ehc) and printed. nag mv cluster indicator finds two clusters and the results are printed. void test_nag_mv_cluster_indicator() { double cd[9], d[45], dord[10], s[10]; matrix x; x.SetSize(10, 10); double dmin_; double dstep, ydist; double dlevel; int ic[10], ilc[9], iord[10], isx[10],iuc[9]; int nsym; int i, j, k; int m = 3, n = 5; int int_method; int tdx = 10; char **c = NULL; char name[10] = "ABCDE"; char char_dist; char char_scale; char char_update; Nag_ClusterMethod method; Nag_MatUpdate update; Nag_DistanceType dist; Nag_VarScaleType scale; int_method = 5; if (int_method == 1) method = Nag_SingleLink; else if (int_method == 2) method = Nag_CompleteLink; else if (int_method == 3) method = Nag_GroupAverage; else if (int_method == 4) method = Nag_Centroid; else if (int_method == 5) method = Nag_Median; else method = Nag_MinVariance; char_update = 'I'; if (char_update == 'U') update = Nag_MatUp; else update = Nag_NoMatUp; char_dist = 'S'; if (char_dist == 'A') dist = Nag_DistAbs; else if (char_dist == 'E') dist = Nag_DistEuclid; else dist = Nag_DistSquared; char_scale = 'U'; if (char_scale == 'S') scale = Nag_VarScaleStd; else if (char_scale == 'R') scale = Nag_VarScaleRange; else if (char_scale == 'G') scale = Nag_VarScaleUser; else scale = Nag_NoVarScale; x[0][0] = 1; x[0][1] = 5.0; x[0][2] = 2.0; x[1][0] = 2; x[1][1] = 1.0; x[1][2] = 1.0; x[2][0] = 3; x[2][1] = 4.0; x[2][2] = 3.0; x[3][0] = 4; x[3][1] = 1.0; x[3][2] = 2.0; x[4][0] = 5; x[4][1] = 5.0; x[4][2] = 0.0; for (i = 0; i < m; ++i) isx[i] = 1; isx[0] = 0; for (i = 0; i < m; ++i) s[i] = 1.0; k = 2; dlevel = 0.0; nag_mv_distance_mat(update, dist, scale, n, m, x, tdx, isx, s, d); nag_mv_hierar_cluster_analysis(method, n, d, ilc, iuc, cd, iord, dord); printf("\nDistance Clusters Joined\n\n"); for (i = 0; i < n-1; ++i) { printf("%10.3f ",cd[i]); printf("%3c",name[ilc[i]-1]); printf("%3c",name[iuc[i]-1]); printf("\n"); } nsym = 20; dmin_ = 0.0; dstep = cd[n - 2] / (1.0* nsym); nag_mv_dendrogram(Nag_DendSouth, n, dord, dmin_, dstep, nsym, &c); printf("\n"); printf("Dendrogram "); printf("\n"); printf("\n"); ydist = cd[n - 2]; for (i = 0; i < nsym; ++i) { if ((i+1) % 3 == 1) { printf("%10.3f%6s",ydist,""); printf("%s",c[i]); printf("\n"); } else { printf("%16s%s","", c[i]); printf("\n"); } ydist -= dstep; } printf("\n"); printf("%14s",""); for (i = 0; i < n; ++i) { printf("%3c",name[iord[i]-1]); printf(""); } printf("\n"); nag_mv_cluster_indicator(n, cd, iord, dord, &k, &dlevel, ic); printf("\n%s%2i%s\n\n","Allocation to ",k," clusters"); printf("Object Cluster\n\n"); for (i = 0; i < n; ++i) { printf("%5s%c%5s","",name[i],""); printf("%ld ",ic[i]); printf("\n"); } } The output is as following: Distance Clusters Joined 1.000 B D 2.000 A C 6.500 A E 14.125 A B Dendrogram 14.125 ------- I I I I 12.006 I I I I I I 9.887 I I I I I I 7.769 I I ---* I I I I 5.650 I I I I I I I I I 3.531 I I I I I I ---* I I 1.412 I I I ---* I I I I I A C E B D Allocation to 2 clusters Object Cluster A 1 B 2 C 1 D 2 E 1 Return: The functions returns NAG error code, 0 if no error. 11: On entry, n must not be less than 2: n = _value_. 19: On entry, k = _value_ while n = _value_. These parameters must satisfy k = n. 118: On entry, dlevel = _value_, k = _value_. Constraint: k = 0 and dlevel > 0.0. 62: The sequence cd is not increasing: cd[_value_] = _value_, cd[_value_] = _value_. 122: On entry, dlevel = _value_, cd[_value_] = _value_. Trivial solution returned. 120: On exit, k = 1 . Trivial solution returned. 121: On exit, k = _value_, n = _value_. Trivial solution returned. 116: Arrays cd and dord are not compatible. 528: The precise number of clusters requested is not possible because of tied clustering distances. The actual number of clusters produced is _value_. 74: An internal error has occurred in this function. Check the function call and array sizes. If the function call is correct, please consult NAG for assistance. successfully call of the nag_mv_cluster_indicator function. */ int nag_mv_cluster_indicator( int n, // the number of objects. double cd[], // the clustering distances in increasing order as returned by nag mv hierar cluster analysis (g03ecc). int iord[], // the objects in the dendrogram order as returned by nag mv hierar cluster analysis double dord[], // the clustering distances corresponding to the order in iord. int *k, // Input: indicates if a speci.ed number of clusters is required. Output: the number of clusters produced, k. double *dlevel, // Input: if k = 0, then dlevel must contain the distance at which clusters are produced. Output: the distance at which the required number of clusters are found if k > 0. int ic[] ); /** g03fac performs a principal coordinate analysis also known as classical metric scaling. Example: The data, given by Krzanowski (1990), are dissimilarities between water vole populations in Europe. The first two principal co-ordinates are computed by nag mv prin coord analysis. void test_nag_mv_prin_coord_analysis() { double e[14]; matrix x; x.SetSize(14, 14); int ndim; int tdx=14; int i, j, n; int nn; char char_roots; Nag_Eigenvalues roots; double d[91] = {0.099, 0.033, 0.022, 0.183, 0.114, 0.042, 0.148, 0.224, 0.059, 0.068, 0.198, 0.039, 0.053, 0.085, 0.051, 0.462, 0.266, 0.322, 0.435, 0.268, 0.025, 0.628, 0.442, 0.444, 0.406, 0.240, 0.129, 0.014, 0.113, 0.070, 0.046, 0.047, 0.034, 0.002, 0.106, 0.129, 0.173, 0.119, 0.162, 0.331, 0.177, 0.039, 0.089, 0.237, 0.071, 0.434, 0.419, 0.339, 0.505, 0.469, 0.390, 0.315, 0.349, 0.151, 0.430, 0.762, 0.633, 0.781, 0.700, 0.758, 0.625, 0.469, 0.618, 0.440, 0.538, 0.607, 0.530, 0.389, 0.482, 0.579, 0.597, 0.498, 0.374, 0.562, 0.247, 0.383, 0.387, 0.084, 0.586, 0.435, 0.550, 0.530, 0.552, 0.509, 0.369, 0.471, 0.234, 0.346, 0.456, 0.090, 0.038}; n = 14; ndim = 2; char_roots = 'L'; nn = n * (n - 1) / 2; if (char_roots == 'L') roots = Nag_LargeEigVals; else roots = Nag_AllEigVals; nag_mv_prin_coord_analysis(roots, n, d, ndim, x, tdx, e); printf("\nScaled Eigenvalues\n\n"); if (char_roots == 'L') { for (i = 0; i < ndim; ++i) printf("%10.4f",e[i]); } else { for (i = 0; i < n; ++i) printf("%10.4f",e[i]); printf("\n"); } printf("\nCo-ordinates\n\n"); for (i = 0; i < n; ++i) { for (j = 0; j < ndim; ++j) printf("%10.4f",x[i][j]); printf("\n"); } } The output is as following: Scaled Eigenvalues 0.7871 0.2808 Co-ordinates 0.2408 0.2337 0.1137 0.1168 0.2394 0.0760 0.2129 0.0605 0.2495 -0.0693 0.1487 -0.0778 -0.0514 -0.1623 0.0115 -0.3446 -0.0039 0.0059 0.0386 -0.0089 -0.0421 -0.0566 -0.5158 0.0291 -0.3180 0.1501 -0.3238 0.0475 Return: The function returns NAG error code or 0 if no error. 70: On entry, parameter roots had an illegal value. 11: On entry, ndim must not be less than 1: ndim = _value_. 18: On entry, n = _value_ while ndim = _value_. These parameters must satisfy n > ndim. 17: On entry, tdx = _value_ while ndim = _value_. These parameters must satisfy tdx = ndim. 115: On entry, d[_value_] = _value_. Constraint: d[i - 1] = 0.0, i = 1, 2, . . . , n * (n - 1)/2. 529: There are fewer than ndim eigenvalues greater than zero. Try a smaller number of dimensions (ndim) or use non-metric scaling (nag mv ordinal multidimscale (g03fcc)). 530: The computation of eigenvalues or eigenvectors has failed. 73: Memory allocation failed. 74: An internal error has occurred in this function. Check the function call and any array sizes. If the call is correct then please consult NAG for assistance. successfully call of the nag_mv_prin_coord_analysis function. */ int nag_mv_prin_coord_analysis( Nag_Eigenvalues roots, int n, // the number of objects in the distance matrix. double d[], // the lower triangle of the distance matrix D stored packed by rows. int ndim, // the number of dimensions used to represent the data. double x[], // the ith row contains k coordinates for the ith point, i = 1, 2, . . . , n. int tdx, // the last dimension of the array x. double eval[] // if roots = Nag AllEigVals, eval contains the n scaled eigenvalues of the matrix E. If roots = Nag LargeEigVals, eval contains the largest k scaled eigenvalues of the matrix E. ); /* g03fcc performs non-metric (ordinal) multidimensional scaling. Example The data, given by Krzanowski (1990), are dissimilarities between water vole populations in Europe. Initial estimates are provided by the first two principal co-ordinates computed by nag mv prin coord analysis (g03fac). The two dimension solution is computed using nag mv ordinal multidimscale. //note: this one requires the struct with pointer. // Has problem on this example code void test_nag_mv_ordinal_multidimscale() { #define NMAX 14 #define MMAX 2 #define NNMAX NMAX*(NMAX-1)/2 #define X(I,J) x[(I-1)*NMAX + (J-1)] #define XTMP(I) xtmp[(I)-1] #define YTMP(I) ytmp[(I)-1] double dfit[364], wk[511]; double x[196]; double stress; int ndim; int i, j, n; int nn; int tdx = NMAX; char char_type; Nag_ScaleCriterion type; double d[91] = {0.099, 0.033, 0.022, 0.183, 0.114, 0.042, 0.148, 0.224, 0.059, 0.068, 0.198, 0.039, 0.053, 0.085, 0.051, 0.462, 0.266, 0.322, 0.435, 0.268, 0.025, 0.628, 0.442, 0.444, 0.406, 0.240, 0.129, 0.014, 0.113, 0.070, 0.046, 0.047, 0.034, 0.002, 0.106, 0.129, 0.173, 0.119, 0.162, 0.331, 0.177, 0.039, 0.089, 0.237, 0.071, 0.434, 0.419, 0.339, 0.505, 0.469, 0.390, 0.315, 0.349, 0.151, 0.430, 0.762, 0.633, 0.781, 0.700, 0.758, 0.625, 0.469, 0.618, 0.440, 0.538, 0.607, 0.530, 0.389, 0.482, 0.579, 0.597, 0.498, 0.374, 0.562, 0.247, 0.383, 0.387, 0.084, 0.586, 0.435, 0.550, 0.530, 0.552, 0.509, 0.369, 0.471, 0.234, 0.346, 0.456, 0.090, 0.038}; n = 14; ndim = 2; char_type = 'T'; nn = n * (n - 1) / 2; nag_mv_prin_coord_analysis(Nag_LargeEigVals, n, d, ndim, x, tdx, wk); if (char_type == 'T') type = Nag_Stress; else type = Nag_SStress; g03fcc(type, n, ndim, d, x, tdx, &stress, dfit, E04_DEFAULT, NAGERR_DEFAULT); Vprintf("\n STRESS = %13.4e\n\n",stress); Vprintf("Co-ordinates\n\n"); for (i = 1; i <= n; ++i) { for (j = 1; j <= ndim; ++j) Vprintf("%10.4f",X(i,j)); 3.g03fcc.4 [NP3275/5/pdf] g03 - Multivariate Methods g03fcc Vprintf("\n"); } exit(EXIT_SUCCESS); } else { Vprintf("Incorrect input value of n.\n"); exit(EXIT_FAILURE); } } } Return: The function returns NAG error code or 0 if no error. 70: On entry, parameter type had an illegal value. 11: On entry, ndim must not be less than 1: ndim = _value_. 18: On entry, n = _value_ while ndim = _value_. These parameters must satisfy n > ndim. 17: On entry, tdx = _value_ while ndim = _value_. These parameters must satisfy tdx = ndim. 531: All elements of array d = 0.0. Constraint: At least one element of d must be positive. 73: Memory allocation failed. 74: An internal error has occurred in this function. Check the function call and array sizes. If the function call is correct, please consult NAG for assistance. Additional error messages are output if the optimization fails to converge or if the options are set incorrectly, Details of these can be found in the nag opt conj grad (e04dgc) document. int nag_mv_ordinal_multidimscale( Nag_ScaleCriterion type, int n, // the number of objects in the distance matrix. int ndim, // the number of dimensions used to represent the data. double d[], // the lower triangle of the distance matrix D stored packed by rows. double x[], // Input: the ith row must contain an initial estimate of the co-ordinates for the ith point, i = 1, 2, . . . , n. One method of computing these is to use nag mv prin coord analysis (g03fac). Output: the ith row contains m co-ordinates for the ith point, i = 1, 2, . . . , n. int tdx, // the last dimension of the array x as declared in the calling program. double *stress, // the value of STRESS or SSTRESS at the final iteration. double dfit[], // auxiliary outputs. Nag_E04_Opt *options, ); */ /** g03xzc frees Nag allocated memory for the dendrogram array in nag mv dendrogram (g03ehc). Example: For an example of the use of nag mv dend free, please refer to the example program for nag mv dendrogram (g03ehc). Return: The function returns NAG error code or 0 if no error. */ int nag_mv_dend_free( char ***c // Input: a pointer to an array of character pointers to which memory has been allocated internally by nag mv dendrogram (g03ehc). Output: The memory allocated to each of the pointers in the pointer array is freed and the pointers are set to NULL. ); /** g03zac produces standardized values (z-scores) for a data matrix. Example: A 4 by 3 data matrix is input along with location and scaling values. The first and third columns are scaled and the results printed. void test_nag_mv_z_scores() { double e[3] = {14.75, 0.0, 1050.0}, s[3] = {2.50, 0.0, 420.0}; matrix z; z.SetSize(4,3); matrix x = {{15.0, 0.0, 1500.0}, {12.0, 1.0, 1000.0}, {18.0, 2.0, 1200.0}, {14.0, 3.0, 500.0}}; int nvar = 2; int isx[3]; int i, j, m = 3, n = 4; int tdx = 3; int tdz = 3; for (j = 0; j < m; ++j) isx[j] = 1; isx[1] = 0; nag_mv_z_scores(n, m, x, tdx, nvar, isx, s, e, z, tdz); printf("\nStandardized Values\n\n"); for (i = 0; i < n; ++i) { for (j = 0; j < nvar; ++j) printf("%8.3f",z[i][j]); printf("\n"); } } The output is following: Standardized Values 0.100 1.071 -1.100 -0.119 1.300 0.357 -0.300 -1.310 Return: The function returns NAG error code or 0 if no error. 11: On entry, n must not be less than 1: n = _value_. On entry, nvar must not be less than 1: nvar = _value_. 17: On entry, m = _value_ while nvar = _value_. These parameters must satisfy m = nvar. On entry, tdx = _value_ while m = _value_. These parameters must satisfy tdx = m. On entry, tdz = _value_ while nvar = _value_. These parameters must satisfy tdz = nvar. 501: The number of variables, nvar in the analysis = _value_, while number of variables included in the analysis via array isx = _value_. Constraint: these two numbers must be the same. 123: On entry, isx[_value_] = _value_, s[_value_] = _value_. Constraint: if isx[j - 1] =0 then s[j - 1] > 0.0, j=1,2,. . .,m. 74: An internal error has occurred in this function. Check the function call and array sizes. If the function call is correct, please consult NAG for assistance. successfully call of the nag_mv_z_scores function. */ int nag_mv_z_scores( int n, // the number of observations in the data matrix. int m, // the number of variables in the data array x. double x[], // contain the ith sample point for the jth variable x[i][j] , for i = 1, 2,. , n, j = 1, 2, . ,m. int tdx, // the last dimension of the array x. int nvar, // the number of variables to be standardised. int isx[], // indicates whether or not the observations on the jth variable are included in the matrix of standardized values. double s[], // the scaling (standard deviation), sj, for the jth variable. double e[], // the location shift (mean), µj, for the jth variable. double z[], // the matrix of standardized values (z-scores). int tdz // the last dimension of the array z ); /* end proto */ #endif /* not NAGG03 */