HW9-solutions

.html

School

University of Texas *

*We aren’t endorsed by this school

Course

230

Subject

Statistics

Date

Apr 3, 2024

Type

html

Pages

Uploaded by sjobs3121

Homework 9 for ISEN 355 (System Simulation) ¶ (C) 2023 David Eckman Due Date: Upload to Canvas by 9:00pm (CDT) on Monday, April 10. In [ ]: # Import some useful Python packages. import numpy as np import matplotlib.pyplot as plt import scipy.stats Problem 1. (10 points) ¶ (a) (2 points) Recall the linear congruential generator (LCG) introduced in Lecture 11. Write a function that takes as input parameters $m$, $a$, and $c$ and a current value $X_n$ and returns the next random number in the sequence, $X_{n+1}$. In Python, the ``%'' symbol acts as the modulus operator. See Section 6.6 of the Python documentation . In [ ]: # Example of modulus with % symbol. # 5 mod 3 = 2 5 % 3 Out[ ]: 2 In [ ]: def myLCG(current_x, m, a, c): next_x = (a * current_x + c) % m return next_x (b) (5 points) As we did in the LCG example in lecture, set $m=32$, $a=11$, and $c=0$. An LCG with $m=32$ has the potential to generate numbers in the set $\{0, 1, 2, \ldots, 30, 31\}$, but as we saw in class, for this choice of $a$ and $c$, the period is less than 32. For a fixed seed, e.g., $X_0 = 1$, recursively run your algorithm until the sequence of returned numbers repeats. Take turns trying different seeds until you have discovered all cycles of numbers. List all cycles. (Note: Each number between 0 and 31 should appear exactly once in your list of cycles.) In [ ]: m = 32 a = 11 c = 0 ListOfCycles = [] repeat = set() for seed in range(m): period = 0 if seed in repeat: continue # If we've seen this number already, skip it. CurrentCycle = [] current_x = seed while current_x not in repeat: repeat.add(current_x) period = period + 1 CurrentCycle.append(current_x) current_x=myLCG(current_x, m, a, c) ListOfCycles.append(CurrentCycle) print(f"The cycle {CurrentCycle} has a period of {period}.")

The cycle [0] has a period of 1. The cycle [1, 11, 25, 19, 17, 27, 9, 3] has a period of 8. The cycle [2, 22, 18, 6] has a period of 4. The cycle [4, 12] has a period of 2. The cycle [5, 23, 29, 31, 21, 7, 13, 15] has a period of 8. The cycle [8, 24] has a period of 2. The cycle [10, 14, 26, 30] has a period of 4. The cycle [16] has a period of 1. The cycle [20, 28] has a period of 2. (c) (3 points) Repeat part (b) with the values $m=32$, $a=13$, and $c=3$. What is the period of the generator for these choices of parameters? In [ ]: m = 32 a = 13 c = 3 ListOfCycles = [] repeat = set() for seed in range(m): period = 0 if seed in repeat: continue # If we've seen this number already, skip it. CurrentCycle = [] current_x = seed while current_x not in repeat: repeat.add(current_x) period = period + 1 CurrentCycle.append(current_x) current_x = myLCG(current_x, m, a, c) ListOfCycles.append(CurrentCycle) print(f"The cycle {CurrentCycle} has a period of {period}.") The cycle [0, 3, 10, 5, 4, 23, 14, 25, 8, 11, 18, 13, 12, 31, 22, 1, 16, 19, 26, 21, 20, 7, 30, 9, 24, 27, 2, 29, 28, 15, 6, 17] has a period of 32. Problem 2. (5 points) ¶ For your simulation model of ETB, how many random numbers do you expect you might need to simulate one day of activity in the building? You may assume that for each random variate needed by your simulation, it can be transformed from one Uniform(0, 1) random variate. Explain in detail how you arrived at your answer. A simulation model of ETB would need to generate random numbers for the interarrival times of students, the schedules of students, the attendance of students to the classes, among other things. Using the data we have, we can obtain the expected number of individuals entering ETB each day; call this quantity X. To generate the arrival times of the X students, we would need at least X random numbers. (If using inversion to generate arrivals from a NSPP, which we didn't discuss in class, this would be true. If using thinning, the number of random numbers needed would depend on the peak of the arrival rate function, because we generate arrivals at the fastest rate.) For these X students, we might need to generate their daily class schedules, which would feature a random number of classes (likely less than or equal to 4) and random selections of classes. As an upper bound, we might need as many as 5X random numbers: 5 for each student (1 to generate a number 1-4 and then that many of numbers to pick the classes). We also need random numbers to handle the choice of which entrance to use ( 1X), taking the elevators or stairs (upper bounded by 1X), and

how long to spend studying (upper bounded by 3X, assuming no one has more than 3 study periods in a day). If the system is more complicated, more variates might be needed, such as for the activities of non-students (e.g., whether they take an elevator, to which floor they go, how long they spend in the building). From the analysis above, we might be looking at around 11X random numbers. If X is around 1500, for example, we might need about 16500 random numbers to simulate one day of operations in ETB. Answers will vary greatly depending on what level of detail you intend to have in your Simio model. Problem 3. (12 points) ¶ (a) (3 points) Print the first 5 Uniform[0, 1] random variates generated by scipy.stats.uniform.rvs when the underlying pseudo-random number generator is seeded with initial state 1234. Read the documentation of the scipy.stats.uniform class and the random_state argument to see how the size and random_state arguments should be set. In [ ]: random_variates = scipy.stats.uniform.rvs(size=5, random_state=1234) print(f"The first 5 random variates are: {random_variates}.") The first 5 random variates are: [0.19151945 0.62210877 0.43772774 0.78535858 0.77997581]. (b) (4 points) With the generator again seeded at state 1234, generate 10,000 Uniform[0, 1] random variates and plot a histogram of their values using 100 bins. Comment on the shape of the histogram. In [ ]: random_variates = scipy.stats.uniform.rvs(size=10000, random_state=1234) plt.hist(random_variates, bins=100) plt.xlabel(r"Random Variate ($U$)") plt.ylabel("Frequency") plt.title("Histogram of 10,000 Uniform[0, 1] RVs") plt.show()

Your preview ends here

Eager to read complete document? Join bartleby learn and gain access to the full version