#################################################################################

# example of transformations of random variables

# unbuffer the output

###########################################

# (1) simulate from a distribution and look at
#     the behavior of transformed values 

# set the number of simulation replications to a large value
# (here i'm using 10 million)

M <- 10000000

# initialize the random number stream for replicability

set.seed( 3773718 )

# generate M IID draws from the PDF of X, namely Uniform( -1, +1 )

x.star <- runif( M, -1, 1 )

# visualize the PDF of X

par( mfrow = c( 1, 2 ) )

hist( x.star, prob = T, breaks = 100, main = 'PDF of X',
  xlab = 'simulated X random draws', xlim = c( -1, 1 ) )

abline( h = 0.5, lwd = 2, col = 'red' )

# transform the simulated X values into simulated Y values

y.star <- x.star^2

# visualize the PDF of Y

hist( y.star, prob = T, breaks = 100, main = 'PDF of Y',
  xlab = 'simulated Y random draws', xlim = c( 0, 1 ) )

y.grid <- seq( 1E-6, 1, length = 500 )

lines( y.grid, 1 / ( 2 * sqrt( y.grid ) ), lwd = 2, col = 'red' )

# let's have a small brush with infinity

hist( y.star, prob = T, breaks = 100, main = 'PDF of Y',
  xlab = 'simulated Y random draws', xlim = c( 0, 1 ),
  ylim = c( 0, 33 ) )

hist( y.star, prob = T, breaks = 1000, main = 'PDF of Y',
  xlab = 'simulated Y random draws', xlim = c( 0, 1 ),
  ylim = c( 0, 33 ) )

# focusing in on the interesting range from 0 to (say) 0.02
# in the right-hand panel

hist( y.star, prob = T, breaks = 1000, main = 'PDF of Y',
  xlab = 'simulated Y random draws', xlim = c( 0, 1 ),
  ylim = c( 0, 1550 ) )

segments( -0.05, 33, 0.05, 33, lwd = 2, col = 'red' )

text( 0.1, 33, '33', col = 'red', cex = 0.75 )

hist( y.star[ y.star < 0.02 ], prob = T, breaks = 1000, main = 'PDF of Y',
  xlab = 'simulated Y random draws', xlim = c( 0, 0.02 ) )

par( mfrow = c( 1, 1 ) )

###########################################

# (2) use the inverse probability integral transform to
#     simulate pseudo-random draws from a target distribution

# start over with an R session with no objects defined

rm( list = ls( ) )

# set a new random number seed for replicability

set.seed( 2385529 )

# set the number of simulation replications to a large value
# (here i'm using 10 million)

M <- 10000000

# generate M IID uniform( 0, 1 ) pseudo-random draws

system.time( 

  u.star <- runif( M )

)

# look at the distribution of the u.star values

hist( u.star, prob = T, breaks = 100, main = 'PDF of U',
  xlab = 'simulated Uniform random draws', xlim = c( 0, 1 ) )

abline( h = 1, lwd = 2, col = 'red' )

# compute the transformed X values

x.star <- u.star / ( 1 - u.star )

# look at the resulting distribution and compare it
# with its theoretical PDF

hist( x.star, prob = T, breaks = 1000, main = 'PDF from voltage example',
  xlab = 'simulated X values' )

# interesting: a few of the u.star values must have come out
# really close to 1

mean( x.star > 25 )

# [1] 0.0384797

# the following plot omits about 3.8% of the distribution
# in the right tail, to concentrate on the interesting part
# of the PDF

hist( x.star[ x.star < 25 ], prob = T, breaks = 1000, 
  main = 'PDF from voltage example', xlab = 'simulated X values' )

x.grid <- seq( 0, 25, length = 500 )

lines( x.grid, 1 / ( 1 + x.grid )^2, lwd = 2, col = 'red' )

#################################################################################