Discussion Problem

Attempt to implement a sampler for 9.3 or 9.4 in Rizzo’s text.

MCMC Methods

Markov Chain Monte Carlo (MCMC) methods encompass a general framework of methods introduced by Metropolis et al. and Hastings for Monte Carlo integration. Many applications of Markov Chain Monte Carlo methods are problems that arise in Bayesian inference because in practice, the normalizing constant for posterior densities are often difficult to evaluate. Markov Chain Monte Carlo provides a method for integrating mathematically intractable equations that are difficult to compute by numerical methods, especially in higher dimensions.

Problem 9.3

Use the Metropolis-Hastings sampler to generate random variables from a standard Cauchy distribution. Discard the first 1000 of the chain, and compare the deciles of the generated observations with the deciles of the standard Cauchy distribution (see qcauchy or qt with df=1). Recall that a Cauchy(\(\theta,\eta\)) distribution has density function \(f(x)= \pi \theta \left[ { 1+\left( { \left[ x-\eta \right] }/{ \theta } \right) }^{ 2 } \right] ^{ -1 }, -\infty >x>\infty ,\theta >0\) The standard Cauchy has the Cauchy(\(\theta=1,\eta=0\)) density. (Note that the standard Cauchy density is equal to the Student \(t\) density with one degree of freedom.)

\[\textrm{Target Distribution:} \quad\quad \frac { 1 }{ \pi \theta \left[ { 1+\left( \frac { x-\eta }{ \theta } \right) }^{ 2 } \right] } ;\theta =1,\eta =0\quad \equiv \quad \frac { \Gamma \left( \frac { v+1 }{ 2 } \right) }{ \sqrt { v\pi } \Gamma \left( \frac { v }{ 2 } \right) } { \left( 1+\frac { { x }^{ 2 } }{ v } \right) }^{ -\frac { v+1 }{ 2 } };v=1\] Generate \(X_{0}\) from the target distribution \(f(x)\).

X <- c(rt(1, df=1))

Generate subsequent states by repeating algorithm in a loop until the chain converges to stationary distribution. The choice of proposal distribution is very flexible, but the chain generated by this choice must satisfy certain regularity conditions. The proposal distribution must be chosen so that the generated chain will converge to a stationary distribution (the target distribution \(f\)). Required conditions for the generated chain are irreducibility, positive recurrence, and aperiodicity

for (t in 2:10^4) {
  Xt <- X[t-1]
  # Generate Y from g(?|X_{t}).
  Y <- rnorm(1)
  # Generate U from Uniform(0,1).
  U <- runif(1)
  num <- dt(Y, df=1) * dnorm(Xt)
  den <- dt(Xt, df=1) * dnorm(Y)
  if (U <= num/den) {
    # accept Y and set X_{t+1} = Y.
    X[t] <- Y 
    }
  else {
    # reject Y and set X_{t+1} = X_{t}.
    X[t] <- Xt
    }
}

X <- X[1001:length(X)]
chain <- quantile(X, seq(0, 1, 0.1))
cauchy <- qcauchy(seq(0, 1, 0.1))
cbind(cauchy, chain)
##          cauchy       chain
## 0%         -Inf -3.93856174
## 10%  -3.0776835 -2.06240245
## 20%  -1.3763819 -1.11420734
## 30%  -0.7265425 -0.65389978
## 40%  -0.3249197 -0.30161255
## 50%   0.0000000 -0.02490065
## 60%   0.3249197  0.24301769
## 70%   0.7265425  0.55568166
## 80%   1.3763819  0.99531437
## 90%   3.0776835  1.69389911
## 100%        Inf  3.42605212

Problem 9.4

Implement a random walk Metropolis sampler for generating the standard Laplace distribution (\(\frac { 1 }{ 2 } { e }^{ \left| x \right| },x\in { R }\)). For the increment, simulate from a normal distribution. Compare the chains generated when different variances are used for the proposal distribution. Also, compute the acceptance rates of each chain.

rw.Metropolis <- function(sigma, X0, N) {
  X <- c(X0)
  k <- 0
  for (i in 2:N) {
      Y <- rnorm(1, X[i-1], sigma)
      U <- runif(1)
      if (U <= (exp(abs(Y)) / exp(abs(X[i-1])))) {
        X[i] <- Y 
        }
      else {
        X[i] <- X[i-1]
        k <- k + 1
        }
    }
  return(list(X=X, k=k))
}

N <- 2000
sigma <- c(.05, .25, .5, .75)
X0 <- 25

rw1 <- rw.Metropolis(sigma[1], X0, N)
rw2 <- rw.Metropolis(sigma[2], X0, N)
rw3 <- rw.Metropolis(sigma[3], X0, N)
rw4 <- rw.Metropolis(sigma[4], X0, N)

Convergence of the random walk Metropolis is often sensitive to the choice of scale parameter. When variance of the increment is too large, most of the candidate points are rejected and the algorithm is very inefficient. If the variance of the increment is too small, the candidate points are almost all accepted, so the random walk Metropolis generates a chain that is almost like a true random walk, which is also inefficient. Below chains are generated with different variances in the proposal distribution.

chain1 <- quantile(rw1$X, seq(0, 1, 0.1))
chain2 <- quantile(rw2$X, seq(0, 1, 0.1))
chain3 <- quantile(rw3$X, seq(0, 1, 0.1))
chain4 <- quantile(rw4$X, seq(0, 1, 0.1))
cbind(chain1, chain2, chain3, chain4)
##        chain1   chain2    chain3    chain4
## 0%   24.23803 24.96403  25.00000  23.76851
## 10%  24.83376 32.87617  47.59344  63.21485
## 20%  25.10013 37.69663  69.06892 102.41206
## 30%  25.32364 45.94662  88.67661 137.04522
## 40%  26.41253 51.90976 122.96196 168.38736
## 50%  26.62386 56.04026 141.41435 208.57444
## 60%  26.89756 62.87835 153.26797 238.72326
## 70%  27.23311 66.09122 175.69053 261.81092
## 80%  27.55548 75.40580 195.70214 288.44252
## 90%  28.04120 83.89074 216.31549 331.54376
## 100% 28.52243 88.43658 233.29124 365.19422

Number of candidate points rejected.

c(rw1$k, rw2$k, rw3$k, rw4$k)
## [1]  31 169 283 393

Assignment Problems

Kelton Chapter 8 problem 5

Problem 8.4.5

Speedee Car Rentals is Open 24 hours and has two classes of customers - 25% are premium customers who are promised a waiting time of less than 5 minutes from when they enter or they will receive a $15 discount on their rental. The others are regular customers who will wait as long as they have to wait. Customers arrive about 2 minutes apart (exponential) and move about 4 meters at 1 mile per hour to either the premium counter (1 agent) or the regular counter (2 agents). All customers have a uniformly distributed service time of from 2 to 7 minutes. Each agent costs $55 per hour (including Overhead) and each customer completed contributes $8 to income (minus any discount they might have received). Model this system without using any links. Run for 10 days per replication. For each scenario, report the utilization of each type of agent and the waiting time of each type of Customer. Both on-screen and in the output, report the number of premium customers who received a discount, the total discounts paid, and the net profit (after discounts and costs) of the system.

Initial Setup

Using the AppointmentArrivals.spfx, FreeSpaceMovement.spfx, and Financials.spfx SimBit examples; created a model with the below characteristics. For the $15 discount for waits over 5 minutes from the Random.Uniform( 2 , 7 ) servers, the probability of waiting over 5 was used: \[\int _{ 5 }^{ 7 }{ \frac { 1 }{ 7-2 } } dx=\int _{ 5 }^{ 7 }{ \frac { 1 }{ 5 } } dx={ \frac { x }{ 5 } \vert }_{ 5 }^{ 7 }=\frac { 7 }{ 5 } -\frac { 5 }{ 5 } =\frac { 2 }{ 5 }\]

  • Two Entity objects with:
    • Name of RegularCustomer and PremiumCustomer.
    • Initial Desired Speed of 1 in Units of Mile per Hour.
    • Initial Travel Mode of Free Space Only.
  • Two Source objects with:
    • Name of RegularSource and PremiumSource.
    • Entity Arrival Logic of Entity Type of CustomerData.CustomerType.
    • Entity Arrival Logic of Interarrival Time of Random.Exponential( 2 ).
    • Table Row Referencing Table Name of CustomerData.
    • Table Row Referencing Row Number of CustomerData.CustomerMix.RandomRow.
  • Two Server objects with:
    • External Input Nodes 4 blocks (to represent 4 meters) awy from Server External Output Nodes.
    • Name of PremiumCounter and RegularCounter.
    • Processing Logic Processing Time of Random.Uniform( 2 , 7 ).
    • Processing Logic Initial Capacity of 2 for Regular and 1 for Premium.
    • Financials Parent Cost Center of RegularCostCenter and PremiumCostCenter.
    • Financials Resource Cost > Usage Cost Rate in USD per Hour of:
      • RegularCounter.Capacity.Initial * 55 for Regular.
      • PremiumCounter.Capacity.Initial * 55 for Premium.
    • Financials Buffer Cost > Input Cost > Cost Per Use for Premium of:
      • 15 * (PremiumCounter.InputBuffer.Contents.TimeWaiting( PremiumCustomer ) > 5.0).
  • One Sink object with:
    • Financials Parent Cost Center of RegularIncomeCenter and PremiumIncomeCenter.
    • Financials Buffer Costs > Input Buffer > Cost Per Use of 8.
  • Add Data Table named CustomerData then:
    • Add Object Reference for Entity with name CustomerType.
    • Marked CustomerType with Set Column As Key.
    • Add Standard Property of Real with name CustomerMix.
    • Added rows with {PremiumCustomer, 0.25} and {RegularCustomer, 0.75}.
  • Add Sequence Table named ServiceTable then:
    • Add Foreign Key with Table Key CustomerData.CustomerType and Name of CustomerType.
    • Add Standard Property of Expression with name ServiceTimes.
    • Added rows with station entity and time for each step and path.
  • Added Experiment Response variable:
    • Named RegularUtil with expression RegularCounter.Capacity.ScheduledUtilization / 100.
    • Named PremiumUtil with expression PremiumCounter.Capacity.ScheduledUtilization / 100.
    • Named RegularTIS with expression RegularCustomer.Population.TimeInSystem.Average.
    • Named PremiumTIS with expression PremiumCustomer.Population.TimeInSystem.Average.
    • Named Discounts with expression PremiumCounter.InputBuffer.Contents.AverageNumberWaiting * (2 / 5).
Model Data Tab

Scenario A

Assume that each agent only handles customers of his designated type.

Results

Scenario B

Relax the above assumption. Consider alternate work rules (who services which customer under which conditions) for this same staff. Evaluate each proposed scenario to see the impact on the metrics.

  • Right-click on Initial Capacity and Set Reference Property to Create New Reference Property.
    • Named RegularCapacity for Regular.
    • Named PremiumCapacity for Premium.
Results

Other Scenarios

Model Other possible solutions to improve the system. You may use links if that helps.

Results

References

http://s1.daumcdn.net/editor/fp/service_nc/pencil/Pencil_chromestore.html

C.P. Robert, G. Casella, Introducing Monte Carlo Methods with R, Use R, DOI 10.1007/978-1-4419-1576-4_6, © Springer Science+Business Media,

Simio and Simulation: Modeling, Analysis, Applications 3d Ed. by W. David Kelton, Jeffrey S. Smith and David T. Sturrock with Simio software.

Discrete-Event Systems Simulation, 5th Edition (2010), by Jerry Banks, John S. Carlson, Barry L. Nelson,and David M. Nicol.