# include # include # include # include # include # include # include # include using namespace std; int main ( int argc, char *argv[] ); double *multinormal_sample ( int m, int n, double a[], double mu[] ); void r8mat_print ( int m, int n, double a[], string title ); void r8mat_print_some ( int m, int n, double a[], int ilo, int jlo, int ihi, int jhi, string title ); void r8mat_write ( string output_filename, int m, int n, double table[] ); double *r8po_fa ( int n, double a[] ); double *r8vec_normal_01_new ( int n ); void r8vec_print ( int n, double a[], string title ); void timestamp ( ); //****************************************************************************80 int main ( int argc, char *argv[] ) //****************************************************************************80 // // Purpose: // // normal_dataset() generates a dataset of multivariate normal random values. // // Usage: // // normal_dataset m n mu a // // where // // * M, the spatial dimension, // * N, the number of points to generate, // * MU, the mean vector. // * A, the MxM variance-covariance matrix. // // creates an M by N multivariate normal random dataset and writes it // to the file "normal_M_N.txt"./ // // Licensing: // // This information is distributed under the MIT license. // // Modified: // // 15 May 2024 // // Author: // // John Burkardt // { double *a; string output_filename; int i; int j; int k; int m; double *mu; ostringstream m_ostring; int n; ostringstream n_ostring; double *x; timestamp ( ); srand48 ( time ( NULL ) ); cout << "\n"; cout << "normal_dataset():\n"; cout << " C++ version\n"; cout << " Generate a multivariate normal random dataset.\n"; cout << "\n"; cout << " The program requests input values from the user:\n"; cout << "\n"; cout << " * M, the spatial dimension,\n"; cout << " * N, the number of points to generate,\n"; cout << " * MU, the mean vector of length M.\n"; cout << " * A, the MxM variance-covariance matrix.\n"; cout << "\n"; cout << " The program generates the data, writes it to the file\n"; cout << "\n"; cout << " normal_M_N.txt\n"; cout << "\n"; cout << " where ""M"" and ""N"" are the numeric values specified\n"; cout << " by the user.\n"; // // Get the spatial dimension. // if ( 1 < argc ) { m = atoi ( argv[1] ); } else { cout << "\n"; cout << " Enter the value of M\n"; cin >> m; } cout << "\n"; cout << " Spatial dimension M = " << m << "\n"; // // Get the number of points. // if ( 2 < argc ) { n = atoi ( argv[2] ); } else { cout << "\n"; cout << " Enter the number of points N\n"; cin >> n; } cout << "\n"; cout << " Number of points N = " << n << "\n"; // // Get the mean. // mu = new double[m]; k = 3; if ( 3 < argc ) { for ( i = 0; i < m; i++ ) { mu[i] = atof ( argv[k] ); k = k + 1; } } else { cout << "\n"; cout << " Enter MU:\n"; for ( i = 0; i < m; i++ ) { cin >> mu[i]; } } r8vec_print ( m, mu, " The mean vector M:" ); // // Get the variance-covariance matrix. // a = new double[m*m]; if ( 4 < argc ) { for ( i = 0; i < m; i++ ) { for ( j = 0; j < m; j++ ) { a[i+j*m] = atof ( argv[k] ); k = k + 1; } } } else { cout << "\n"; cout << " Enter A:\n"; for ( i = 0; i < m; i++ ) { for ( j = 0; j < m; j++ ) { cin >> a[i+j*m]; } } } r8mat_print ( m, m, a, " The variance-covariance matrix A:" ); // // Compute the data. // x = multinormal_sample ( m, n, a, mu ); // // Write it to a file. // m_ostring << m; n_ostring << n; output_filename = "normal_" + m_ostring.str ( ) + "_" + n_ostring.str ( ) + ".txt"; r8mat_write ( output_filename, m, n, x ); cout << "\n"; cout << " The data was written to the file \"" << output_filename << "\".\n"; // // Terminate. // delete [] a; delete [] mu; delete [] x; cout << "\n"; cout << "normal_dataset():\n"; cout << " Normal end of execution.\n"; cout << "\n"; timestamp ( ); return 0; } //****************************************************************************80 double *multinormal_sample ( int m, int n, double a[], double mu[] ) //****************************************************************************80 // // Purpose: // // multinormal_sample() samples a multivariate normal distribution. // // Discussion: // // The multivariate normal distribution for the M dimensional vector X // has the form: // // pdf(X) = (2*pi*det(A))**(-M/2) * exp(-0.5*(X-MU)'*inverse(A)*(X-MU)) // // where MU is the mean vector, and A is a positive definite symmetric // matrix called the variance-covariance matrix. // // Licensing: // // This information is distributed under the MIT license. // // Modified: // // 15 May 2024 // // Author: // // John Burkardt // // Input: // // int M, the dimension of the space. // // int N, the number of points. // // double A[M*M], the variance-covariance // matrix. A must be symmetric positive definite. // // double MU[M], the mean vector. // // Output: // // double MULTINORMAL_SAMPLE[M*N], the points. // { int i; int j; int k; double *r; double *x; double *y; // // Compute the upper triangular Cholesky factor R of the variance-covariance // matrix. // r = r8po_fa ( m, a ); if ( !r ) { cout << "\n"; cout << "multinormal_sample(): Fatal error!\n"; cout << " The variance-covariance matrix is not positive definite symmetric.\n"; exit ( 1 ); } // // Y = MxN matrix of samples of the 1D normal distribution with mean 0 // and variance 1. // y = r8vec_normal_01_new ( m*n ); // // Compute X = MU + R' * Y. // x = new double[m*n]; for ( j = 0; j < n; j++ ) { for ( i = 0; i < m; i++ ) { x[i+j*m] = mu[i]; for ( k = 0; k < m; k++ ) { x[i+j*m] = x[i+j*m] + r[k+i*m] * y[k+j*m]; } } } delete [] r; delete [] y; return x; } //****************************************************************************80 void r8mat_print ( int m, int n, double a[], string title ) //****************************************************************************80 // // Purpose: // // r8mat_print() prints an R8MAT. // // Discussion: // // An R8MAT is a doubly dimensioned array of R8 values, stored as a vector // in column-major order. // // Entry A(I,J) is stored as A[I+J*M] // // Licensing: // // This information is distributed under the MIT license. // // Modified: // // 15 May 2024 // // Author: // // John Burkardt // // Input: // // int M, the number of rows in A. // // int N, the number of columns in A. // // double A[M*N], the M by N matrix. // // string TITLE, a title. // { r8mat_print_some ( m, n, a, 1, 1, m, n, title ); return; } //****************************************************************************80 void r8mat_print_some ( int m, int n, double a[], int ilo, int jlo, int ihi, int jhi, string title ) //****************************************************************************80 // // Purpose: // // r8mat_print_some() prints some of an R8MAT. // // Discussion: // // An R8MAT is a doubly dimensioned array of R8 values, stored as a vector // in column-major order. // // Licensing: // // This information is distributed under the MIT license. // // Modified: // // 15 May 2024 // // Author: // // John Burkardt // // Input: // // int M, the number of rows of the matrix. // M must be positive. // // int N, the number of columns of the matrix. // N must be positive. // // double A[M*N], the matrix. // // int ILO, JLO, IHI, JHI, designate the first row and // column, and the last row and column to be printed. // // string TITLE, a title. // { # define INCX 5 int i; int i2hi; int i2lo; int j; int j2hi; int j2lo; cout << "\n"; cout << title << "\n"; // // Print the columns of the matrix, in strips of 5. // for ( j2lo = jlo; j2lo <= jhi; j2lo = j2lo + INCX ) { j2hi = j2lo + INCX - 1; j2hi = min ( j2hi, n ); j2hi = min ( j2hi, jhi ); cout << "\n"; // // For each column J in the current range... // // Write the header. // cout << " Col: "; for ( j = j2lo; j <= j2hi; j++ ) { cout << setw(7) << j << " "; } cout << "\n"; cout << " Row\n"; cout << "\n"; // // Determine the range of the rows in this strip. // i2lo = max ( ilo, 1 ); i2hi = min ( ihi, m ); for ( i = i2lo; i <= i2hi; i++ ) { // // Print out (up to) 5 entries in row I, that lie in the current strip. // cout << setw(5) << i << " "; for ( j = j2lo; j <= j2hi; j++ ) { cout << setw(12) << a[i-1+(j-1)*m] << " "; } cout << "\n"; } } return; # undef INCX } //****************************************************************************80 void r8mat_write ( string output_filename, int m, int n, double table[] ) //****************************************************************************80 // // Purpose: // // r8mat_write() writes an R8MAT file with no header. // // Licensing: // // This information is distributed under the MIT license. // // Modified: // // 15 May 2024 // // Author: // // John Burkardt // // Input: // // string OUTPUT_FILENAME, the output filename. // // int M, the spatial dimension. // // int N, the number of points. // // double TABLE[M*N], the table data. // { int i; int j; ofstream output; // // Open the file. // output.open ( output_filename.c_str ( ) ); if ( !output ) { cerr << "\n"; cerr << "r8mat_write(): Fatal error!\n"; cerr << " Could not open the output file.\n"; return; } // // Write the data. // for ( j = 0; j < n; j++ ) { for ( i = 0; i < m; i++ ) { output << " " << setw(24) << setprecision(16) << table[i+j*m]; } output << "\n"; } // // Close the file. // output.close ( ); return; } //****************************************************************************80 double *r8po_fa ( int n, double a[] ) //****************************************************************************80 // // Purpose: // // r8po_fa() factors a R8PO matrix. // // Discussion: // // The R8PO storage format is appropriate for a symmetric positive definite // matrix and its inverse. (The Cholesky factor of a R8PO matrix is an // upper triangular matrix, so it will be in R8GE storage format.) // // Only the diagonal and upper triangle of the square array are used. // This same storage format is used when the matrix is factored by // R8PO_FA, or inverted by R8PO_INVERSE. For clarity, the lower triangle // is set to zero. // // The positive definite symmetric matrix A has a Cholesky factorization // of the form: // // A = R' * R // // where R is an upper triangular matrix with positive elements on // its diagonal. // // Licensing: // // This information is distributed under the MIT license. // // Modified: // // 15 May 2024 // // Author: // // Original FORTRAN77 version by Dongarra, Bunch, Moler, Stewart. // This version by John Burkardt. // // Reference: // // Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart, // LINPACK User's Guide, // SIAM, 1979, // ISBN13: 978-0-898711-72-1, // LC: QA214.L56. // // Input: // // int N, the order of the matrix. // // double A[N*N], the matrix in R8PO storage. // // Output: // // double R8PO_FA[N*N], the Cholesky factor in SGE // storage, or NULL if there was an error. // { double *b; int i; int j; int k; double s; b = new double[n*n]; for ( j = 0; j < n; j++ ) { for ( i = 0; i < n; i++ ) { b[i+j*n] = a[i+j*n]; } } for ( j = 0; j < n; j++ ) { for ( k = 0; k <= j-1; k++ ) { for ( i = 0; i <= k-1; i++ ) { b[k+j*n] = b[k+j*n] - b[i+k*n] * b[i+j*n]; } b[k+j*n] = b[k+j*n] / b[k+k*n]; } s = b[j+j*n]; for ( i = 0; i <= j-1; i++ ) { s = s - b[i+j*n] * b[i+j*n]; } if ( s <= 0.0 ) { delete [] b; return NULL; } b[j+j*n] = sqrt ( s ); } // // Since the Cholesky factor is in R8GE format, zero out the lower triangle. // for ( i = 0; i < n; i++ ) { for ( j = 0; j < i; j++ ) { b[i+j*n] = 0.0; } } return b; } //****************************************************************************80 double *r8vec_normal_01_new ( int n ) //****************************************************************************80 // // Purpose: // // r8vec_normal_01_new() returns a unit pseudonormal R8VEC. // // Discussion: // // The standard normal probability distribution function (PDF) has // mean 0 and standard deviation 1. // // This routine can generate a vector of values on one call. // // Licensing: // // This information is distributed under the MIT license. // // Modified: // // 09 December 2024 // // Author: // // John Burkardt // // Input: // // int N, the number of values desired. // // Output: // // double R8VEC_NORMAL_01_NEW[N], a sample of the standard normal PDF. // { int i; double r1; double r2; static double r8_pi = 3.141592653589793; double *x; x = new double[n]; for ( i = 0; i < n; i++ ) { r1 = drand48 ( ); r2 = drand48 ( ); x[i] = sqrt ( - 2.0 * log ( r1 ) ) * cos ( 2.0 * r8_pi * r2 ); } return x; } //****************************************************************************80 void r8vec_print ( int n, double a[], string title ) //****************************************************************************80 // // Purpose: // // r8vec_print() prints an R8VEC. // // Discussion: // // An R8VEC is a vector of R8's. // // Licensing: // // This information is distributed under the MIT license. // // Modified: // // 15 May 2024 // // Author: // // John Burkardt // // Input: // // int N, the number of components of the vector. // // double A[N], the vector to be printed. // // string TITLE, a title. // { int i; cout << "\n"; cout << title << "\n"; cout << "\n"; for ( i = 0; i < n; i++ ) { cout << " " << setw(8) << i << " " << setw(14) << a[i] << "\n"; } return; } //****************************************************************************80 void timestamp ( ) //****************************************************************************80 // // Purpose: // // timestamp() prints the current YMDHMS date as a time stamp. // // Example: // // 31 May 2001 09:45:54 AM // // Licensing: // // This information is distributed under the MIT license. // // Modified: // // 15 May 2024 // // Author: // // John Burkardt // { # define TIME_SIZE 40 static char time_buffer[TIME_SIZE]; const struct std::tm *tm_ptr; std::time_t now; now = std::time ( NULL ); tm_ptr = std::localtime ( &now ); std::strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm_ptr ); std::cout << time_buffer << "\n"; return; # undef TIME_SIZE }