polya.sim1 <- function( M, n.sim, cc ) {

  b <- matrix( 0, M, 2^M - 1 )

  for ( i in 1:M ) {

    b[i,1:( 2^i - 1 )] <- qnorm( ( 1:( 2^i - 1 ) ) / 2^i )

  }

  alpha <- matrix( 0, M, 2^M )

  for ( i in 1:M ) {

    alpha[i, 1:( 2^i )] <- cc * rep( i^2, 2^i )

  }

  par( mfrow = c( 2, 1 ) )

  plot( seq( -3, 3, length = 500 ), dnorm( seq( -3, 3, 
    length = 500 ) ), 
    type = 'l', xlab = 'y', ylab = 'Density', ylim = c( 0, 1.5 ) )
#     main = paste( 'c =', cc, ', n.sim =', n.sim ) )

  F.star.cumulative <- rep( 0, 2^M )

  b.star <- c( -3, b[M,], 3 )

  for ( i in 1:n.sim ) {

    F.star <- rep( 1, 2^M )

    for ( j in 1:M ) {

      for ( k in 1:( 2^( j - 1 ) ) ) {

        C <- rbeta( 1, alpha[j, 2 * k - 1], alpha[j, 2 * k] )

        F.star[( 1 + ( k - 1 ) * 2^( M - j + 1 ) ):( ( 2 * k - 1 ) * 
          2^( M - j ) )] <- F.star[( 1 + ( k - 1 ) * 2^( M - j + 1 ) ):
          ( ( 2 * k - 1 ) * 2^( M - j ) )] * C

        F.star[( ( 2 * k - 1 ) * 2^( M - j ) + 1 ):( k * 2^( M - j + 
          1 ) )] <- F.star[( ( 2 * k - 1 ) * 2^( M - j ) + 1 ):( k * 
          2^( M - j + 1 ) )] * ( 1 - C )

      }

    }

    F.star.cumulative <- F.star.cumulative + F.star

    n <- round( 10000 * F.star )

    y <- NULL

    for ( j in 1:2^M ) {

      y <- c( y, runif( n[j], b.star[j], b.star[j+1] ) )

    }

    lines( density( y ), lty = 2 )

    print( i )

  }

  F.star.cumulative <- F.star.cumulative / n.sim

  n <- round( 10000 * F.star.cumulative )

  y <- NULL

  for ( i in 1:2^M ) {

    y <- c( y, runif( n[i], b.star[i], b.star[i+1] ) )

  }

  hist( y, nclass = 20, probability = T, ylab = 'Density',
    xlim = c( -3, 3 ), ylim = c( 0, 0.5 ), xlab = 'y' )

  lines( seq( -3, 3, length = 500 ), dnorm( seq( -3, 3, length = 500 ) ) )

  lines( density( y ), lty = 2 )

  return( cat( "\007" ) )

}

y <- exp( rnorm( 100 ) ) - 2

polya.update <- function( M, n.sim, cc, y ) {

  b <- matrix( 0, M, 2^M - 1 )

  for ( i in 1:M ) {

    b[i,1:( 2^i - 1 )] <- qnorm( ( 1:( 2^i - 1 ) ) / 2^i )

  }

  b.star <- c( -3, b[ M, ], 3 )

  infinity <- 2 * abs( max( min( y ), max( y ), min( b ), max( b ) ) )

  alpha <- matrix( 0, M, 2^M )

  for ( i in 1:M ) {

    alpha[i, 1:( 2^i )] <- cc * rep( i^2, 2^i )

  }

  par( mfrow = c( 1, 1 ) )

  hist( y, xlim = c( -3, max( y ) ), ylim = c( 0, 1.5 ), xlab = 'y', 
    ylab = 'Density', probability = T, nclass = 20,
    main = paste( 'n.sim =', n.sim, ', c =', cc, ', n = ', length( y ) ) )

  for ( i in 1:n.sim ) {

    F.star <- rep( 1, 2^M )

    for ( j in 1:M ) {

      for ( k in 1:( 2^( j - 1 ) ) ) {

        C <- rbeta( 1, alpha[j, 2 * k - 1], alpha[j, 2 * k] )

        F.star[( 1 + ( k - 1 ) * 2^( M - j + 1 ) ):( ( 2 * k - 1 ) * 
          2^( M - j ) )] <- F.star[( 1 + ( k - 1 ) * 2^( M - j + 1 ) ):
          ( ( 2 * k - 1 ) * 2^( M - j ) )] * C

        F.star[( ( 2 * k - 1 ) * 2^( M - j ) + 1 ):( k * 2^( M - j + 
          1 ) )] <- F.star[( ( 2 * k - 1 ) * 2^( M - j ) + 1 ):( k * 
          2^( M - j + 1 ) )] * ( 1 - C )

      }

    }

    n <- round( 10000 * F.star )

    y.star <- NULL

    for ( j in 1:2^M ) {

      y.star <- c( y.star, runif( n[j], b.star[j], b.star[j+1] ) )

    }

    lines( density( y.star ), lty = 1 )

  }

  for ( i in 1:M ) {

    n <- hist( y, breaks = c( - infinity, b[i,1:( 2^i - 1 )], 
      infinity ), plot = F )$counts
    alpha[i, 1:( 2^i )] <- alpha[i, 1:( 2^i )] + n

  }

  for ( i in 1:n.sim ) {

    F.star <- rep( 1, 2^M )

    for ( j in 1:M ) {

      for ( k in 1:( 2^( j - 1 ) ) ) {

        C <- rbeta( 1, alpha[j, 2 * k - 1], alpha[j, 2 * k] )

        F.star[( 1 + ( k - 1 ) * 2^( M - j + 1 ) ):( ( 2 * k - 1 ) * 
          2^( M - j ) )] <- F.star[( 1 + ( k - 1 ) * 2^( M - j + 1 ) ):
          ( ( 2 * k - 1 ) * 2^( M - j ) )] * C

        F.star[( ( 2 * k - 1 ) * 2^( M - j ) + 1 ):( k * 2^( M - j + 
          1 ) )] <- F.star[( ( 2 * k - 1 ) * 2^( M - j ) + 1 ):( k * 
          2^( M - j + 1 ) )] * ( 1 - C )

      }

    }

    n <- round( 10000 * F.star )

    y.star <- NULL

    for ( j in 1:2^M ) {

      y.star <- c( y.star, runif( n[j], b.star[j], b.star[j+1] ) )

    }

    lines( density( y.star ), lty = 2 )

  }

  return( cat( "\007" ) )

}

polya.sim( 8, 100, 10 )

polya.sim( 8, 100, 1 )

polya.sim( 8, 100, 0.1 )

polya.update( 8, 10, 10, y )

polya.update( 8, 10, 1, y )

polya.update( 8, 10, 0.1, y )