par( mfrow = c( 1, 2 ) ) x.grid <- seq( -10, 20, length = 500 ) plot( x.grid, ( x.grid + 10 )^3 / 27000, type = 'l', xlab = 'x', ylab = 'CDF of X', lwd = 2, col = 'red' ) p.grid <- seq( 0.0, 1.0, length = 500 ) plot( p.grid, 30 * p.grid^( 1 / 3 ) - 10, type = 'l', xlab = 'p', ylab = 'Inverse CDF of X', lwd = 2, col = 'red' ) par( mfrow = c( 1, 1 ) ) m <- 100000 set.seed( 8787986 ) u.star <- runif( m ) x.star <- 30 * u.star^( 1 / 3 ) - 10 hist( x.star, breaks = 100, prob = T, main = '', xlab = 'x', ylab = 'simulated PDF of X' ) lines( x.grid, ( x.grid + 10 )^2 / 9000, lty = 1, lwd = 2, col = 'red' )