# this function simulates and plots n.sim draws from a # Polya tree (PT) prior centered at the standard normal # distribution with prior sample size c = cc (also often # known as alpha) and with the theoretically infinite number # of levels of the tree set to the finite value M # (unbuffer the output) set.sim( 2718282 ) 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 ) ) 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" ) ) } # this code plots draws from the PT prior centered at # the standard normal distribution with M = 6, cc = 0.1 and # n.sim = 100 polya.sim1( 6, 100, 0.1 ) # the top panel shows the draws from the PT prior on the # density scale; the bottom panel shows the average of the # PT draws expressed as a histogram and a density trace # (dotted line) and the target centering distribution # (solid line); this demonstrates that, even though the # top panel looks completely chaotic, all of the draws are # actually centered at N( 0, 1 ) (with a larger value of # n.sim the dotted curve would match the solid curve exactly) # now try the same thing with cc = 1 polya.sim1( 6, 100, 1 ) # the draws are visibly closer to the centering distribution # how about cc = 10? polya.sim1( 6, 100, 10 ) # now every draw is close to N( 0, 1 ) # the following code (a) simulates n = 100 data values from # a distribution with a long right-hand tail and (b) samples # and plots draws from the PT posterior that starts with # the PT prior in polya.sim1 and updates with the data in (a) # (y.lim is a hack that scales the plots well) y <- c( rep( -2.75, 24 ), rep( -2.25, 27 ), rep( -1.75, 14 ), rep( -1.25, 13 ), rep( -0.75, 7 ), rep( -0.25, 2 ), rep( 0.25, 2 ), rep( 0.75, 3 ), rep( 1.25, 2 ), rep( 1.75, 1 ), rep( 2.25, 2 ), rep( 2.75, 0 ), rep( 3.25, 1 ), rep( 3.75, 2 ) ) polya.update <- function( M, n.sim, cc, y, y.lim ) { 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( -4, max( y ) ), ylim = c( 0, y.lim ), xlab = 'y', ylab = 'Density', probability = T, nclass = 20, col = 'red', main = paste( 'n.sim =', n.sim, ', c =', cc, ', n = ', length( y ) ) ) standard.normal.grid <- seq( -3, 3, length = 500 ) lines( standard.normal.grid, dnorm( standard.normal.grid ), lwd = 6, col = 'blue' ) 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( "quack!" ) } # let's try running this function with M = 6, n.sim = 100, # cc = 0.1, and the data set y generated above polya.update( 6, 25, 0.1, y, y.lim = 5 ) # the prior draws are in solid black, the data histogram is in red, # the prior centering distribution is in blue, and the posterior # draws are in dotted black # you can see that with c = 0.1 the posterior is coming # almost entirely from the data # same thing with c = 1 polya.update( 6, 25, 1, y, y.lim = 1.5 ) # the posterior is now a compromise between the prior and # the data # same thing with c = 10 polya.update( 6, 25, 10, y, y.lim = 0.6 ) # the posterior is now coming almost entirely from the prior