y <- c( 495, 541, 1461, 1555, 1603, 2201, 2750,
  3468, 3516, 4319, 6622, 7728, 13159, 21194 )

lambda.true <- mean( y )

n <- length( y )

p.grid <- ( ( 1:n ) - 0.5 ) / n

set.seed( 1 )

y.star <- rexp( n, 1 / lambda.true )

lambda.hat.star <- mean( y.star )

plot( - lambda.hat.star * log( 1 - p.grid ), sort( y.star ), type = 'l',
  xlab = 'Exponential Quantiles', ylab = 'Sorted y Values',
  xlim = c( 0, 57500 ), ylim = c( 0, 57500 ), lty = 2 )

M <- 1000

for ( i in 2:M ) {

  y.star <- rexp( n, 1 / lambda.true )

  lambda.hat.star <- mean( y.star )

  lines( - lambda.hat.star * log( 1 - p.grid ), sort( y.star ), lty = 2 )

}

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

lines( - lambda.true * log( 1 - p.grid ), sort( y ), lty = 1,
  lwd = 4, col = 'blue' )