function discrete_pdf_sample_2d_test01 ( n )

%*****************************************************************************80
%
%% discrete_pdf_sample_2d_test01() looks at a 20x20 region.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    15 December 2012
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer N, the number of sample points to be generated.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'discrete_pdf_sample_2d_test01():\n' );
  fprintf ( 1, '  Consider data skewed toward the upper left corner of the unit square.\n' );
  fprintf ( 1, '  Generate %d samples\n', n );
%
%  Get the dimensions of the PDF data.
%
  [ n1, n2 ] = get_discrete_pdf_size1 ( );
  fprintf ( 1, '  PDF data is on a %d by %d grid.', n1, n2 );
%
%  Construct a PDF from the data.
%
  pdf = get_discrete_pdf_data1 ( n1, n2 );
%
%  "Integrate" the data over rows and columns of the region to get the CDF.
%
  cdf = set_discrete_cdf ( n1, n2, pdf );
%
%  Choose N CDF values at random.
%
  u = rand ( n, 1 );
%
%  Find the cell corresponding to each CDF value,
%  and choose a random point in that cell.
%
  xy = discrete_cdf_to_xy ( n1, n2, cdf, n, u );
%
%  Write data to a file for examination, plotting, or analysis.
%
  filename = 'discrete_pdf_sample_2d_test01.txt';
  r8mat_write ( filename, 2, n, xy );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Wrote sample data to file "%s".\n', filename );

  return
end
function pdf = get_discrete_pdf_data1 ( n1, n2 )

%*****************************************************************************80
%
%% get_discrete_pdf_data1() returns the value of the discrete PDF function in each cell.
%
%  Discussion:
%
%    Cell (I,J) extends from 
%
%      (I-1) * H < Y < I * H
%      (J-1) * H < X < J * H
%
%    We have data for each cell, representing the integral of some PDF
%    over that cell.  The function pdf(x,y) must be nonnegative.  However,
%    we don't impose any other conditions on it.
%
%    The array PDF(:,:) contains the integral of pdf(x,y) over each cell,
%    or, almost as good, simply a sample or average value.
%
%    We load the array PDF, and then we normalize it so that the sum of
%    all the entries is 1.  
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    15 December 2012
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer N1, N2, the number of rows and columns of PDF data.
%
%  Output:
%
%    real PDF(N1,N2).  PDF(I,J) is the discrete PDF 
%    for the cell (I,J), normalized so that the sum over all cells is 1.
%

%
%  This is the raw data.
%
  pdf = [ ...
    0.0000,0.0000,0.0000,0.0000,0.0000, ...
    0.0000,0.0000,0.0000,0.0000,0.0000, ...
    0.0000,0.0000,0.0000,0.0000,0.0000, ...
    0.0000,0.0000,0.0000,0.0000,0.0000; ...
    0.0000,0.0000,0.0001,0.0001,0.0002, ...
    0.0002,0.0002,0.0003,0.0003,0.0003, ...
    0.0003,0.0003,0.0002,0.0002,0.0002, ...
    0.0002,0.0001,0.0001,0.0001,0.0000; ...
    0.0000,0.0001,0.0002,0.0003,0.0004, ...
    0.0004,0.0005,0.0006,0.0006,0.0006, ...
    0.0006,0.0006,0.0005,0.0005,0.0004, ...
    0.0003,0.0003,0.0002,0.0001,0.0000; ...
    0.0000,0.0002,0.0003,0.0005,0.0006, ...
    0.0008,0.0009,0.0009,0.0010,0.0010, ...
    0.0010,0.0009,0.0008,0.0008,0.0007, ...
    0.0006,0.0004,0.0003,0.0002,0.0000; ...
    0.0000,0.0003,0.0005,0.0008,0.0010, ...
    0.0012,0.0014,0.0015,0.0015,0.0015, ...
    0.0015,0.0014,0.0013,0.0011,0.0010, ...
    0.0008,0.0006,0.0005,0.0003,0.0000; ...
    0.0000,0.0004,0.0009,0.0013,0.0016, ...
    0.0019,0.0021,0.0023,0.0023,0.0023, ...
    0.0021,0.0020,0.0018,0.0016,0.0013, ...
    0.0011,0.0009,0.0007,0.0004,0.0000; ...
    0.0000,0.0007,0.0014,0.0020,0.0025, ...
    0.0030,0.0033,0.0034,0.0034,0.0033, ...
    0.0031,0.0028,0.0025,0.0022,0.0018, ...
    0.0015,0.0012,0.0009,0.0006,0.0000; ...
    0.0000,0.0011,0.0021,0.0031,0.0039, ...
    0.0045,0.0049,0.0051,0.0050,0.0047, ...
    0.0043,0.0039,0.0034,0.0029,0.0024, ...
    0.0019,0.0015,0.0011,0.0007,0.0000; ...
    0.0000,0.0017,0.0033,0.0048,0.0060, ...
    0.0069,0.0074,0.0074,0.0072,0.0066, ...
    0.0059,0.0052,0.0045,0.0037,0.0031, ...
    0.0025,0.0019,0.0014,0.0009,0.0000; ...
    0.0000,0.0025,0.0050,0.0073,0.0091, ...
    0.0104,0.0109,0.0107,0.0101,0.0091, ...
    0.0080,0.0068,0.0057,0.0047,0.0038, ...
    0.0030,0.0023,0.0017,0.0011,0.0000; ...
    0.0000,0.0038,0.0075,0.0110,0.0136, ...
    0.0153,0.0157,0.0151,0.0138,0.0121, ...
    0.0104,0.0087,0.0071,0.0058,0.0046, ...
    0.0036,0.0027,0.0019,0.0012,0.0000; ...
    0.0000,0.0055,0.0110,0.0160,0.0198, ...
    0.0218,0.0219,0.0205,0.0182,0.0155, ...
    0.0129,0.0106,0.0085,0.0068,0.0053, ...
    0.0041,0.0031,0.0022,0.0014,0.0000; ...
    0.0000,0.0077,0.0154,0.0224,0.0276, ...
    0.0299,0.0293,0.0266,0.0229,0.0190, ...
    0.0154,0.0123,0.0098,0.0077,0.0059, ...
    0.0045,0.0034,0.0024,0.0015,0.0000; ...
    0.0000,0.0100,0.0202,0.0295,0.0362, ...
    0.0385,0.0368,0.0324,0.0271,0.0219, ...
    0.0174,0.0137,0.0107,0.0082,0.0063, ...
    0.0048,0.0035,0.0025,0.0016,0.0000; ...
    0.0000,0.0120,0.0244,0.0356,0.0432, ...
    0.0455,0.0426,0.0366,0.0298,0.0236, ...
    0.0184,0.0143,0.0110,0.0084,0.0064, ...
    0.0048,0.0035,0.0025,0.0016,0.0000; ...
    0.0000,0.0134,0.0266,0.0382,0.0461, ...
    0.0480,0.0445,0.0376,0.0301,0.0235, ...
    0.0181,0.0139,0.0106,0.0081,0.0061, ...
    0.0046,0.0033,0.0023,0.0015,0.0000; ...
    0.0000,0.0151,0.0261,0.0362,0.0436, ...
    0.0447,0.0412,0.0347,0.0276,0.0214, ...
    0.0164,0.0125,0.0095,0.0072,0.0054, ...
    0.0041,0.0029,0.0021,0.0013,0.0000; ...
    0.0000,0.0174,0.0220,0.0295,0.0349, ...
    0.0361,0.0333,0.0281,0.0225,0.0175, ...
    0.0134,0.0102,0.0078,0.0059,0.0044, ...
    0.0033,0.0024,0.0017,0.0010,0.0000; ...
    0.0000,0.0097,0.0152,0.0200,0.0235, ...
    0.0244,0.0227,0.0193,0.0156,0.0122, ...
    0.0094,0.0072,0.0055,0.0041,0.0031, ...
    0.0023,0.0017,0.0012,0.0007,0.0000; ...
    0.0000,0.0000,0.0000,0.0000,0.0000, ...
    0.0000,0.0000,0.0000,0.0000,0.0000, ...
    0.0000,0.0000,0.0000,0.0000,0.0000, ...
    0.0000,0.0000,0.0000,0.0000,0.0000 ]';
%
%  Normalize to get an integral of 1.
%
  total = sum ( sum ( pdf(1:n1,1:n2) ) );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  PDF data sums to %g\n', total );

  pdf(1:n1,1:n2) = pdf(1:n1,1:n2) / total;

  return
end
function [ n1, n2 ] = get_discrete_pdf_size1 ( )

%*****************************************************************************80
%
%% get_discrete_pdf_size1() returns the dimension of the discrete PDF data.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    15 December 2012
%
%  Author:
%
%    John Burkardt
%
%  Output:
%
%    integer N1, N2, the number of rows and columns of data.
%
  n1 = 20;
  n2 = 20;

  return
end
