For this assignment, implement the Bagel Shop problem in Python. Since this is a Monte Carlo simulation (not a Discrete Event Simulation), you won't need to use SimPy.

Using probability concepts, it can be shown that: (1) the expected number of customers per day is 10.2; (2) the expected number of bagels is 2 dozen per customer; and (3) the expected number of bagels needed per day (multiply the answers from (1) and (2)) is 20.4. Based on those values, the owner would make 20 dozen bagels each day, with (4) expected cost of 20*$5.80 = $116 per day, (5) expected income of 20*$8.40 = 168, and (6) an expected profit of $168-$116 = $52.

Use Python to implement a simulation of at least 5 days of operation (call that one week), and perform 400 replications. Investigate the impact of making between 15 and 35 dozen bagels per day on the net weekly profit. How much should the baker make each day?

Write a report that (1) explains your rationale for your recommendation of how many dozens of bagels should be made per day; (2) compares the answer from the simulation to the expected-value solution (above) of 20 dozen; and (3) compares your Python-based simulation results to those from your Excel-based simulation.

News Dealer's Problem from BCNN5, Example 2.7 (pp. 57-60)

Scenario: A news dealer buys papers for $0.33 and sells them for $0.50 each. Unsold newspapers are sold for scrap ($0.05 each). Newspapers are bought in bundles of 10. The demand for newspapers depends on the type of newsday ("good", "fair", or "poor") which is unknown prior to buying the papers. Profit = Revenue - Cost - Lost_profit_from_unmet_demand + Salvage The goal is to determine the order quantity, Q (how many papers to order at the beginning of the newsday).

Created on Sun Jan 31 22:26:00 2016

@author: ThomasR """

#import random from scipy.stats import rv_discrete import numpy as np, scipy.stats as st #import simpy import matplotlib.pyplot as plt #import time.perf_counter as clock_it

class G: #---- Global Variables --- """ To reference a one of these "global" variables, prefix the variable name with a "G." For example: G.REPS or G.PRICE """ REPS = 400 # Number of times to run the monthly sim... PRICE = 0.50 # Selling price COST = 0.33 # Cost to dealer SALVAGE = 0.05 # Value of left over papers ILLWILL = PRICE-COST # lost profit from excess demand DAY_LIST = ["Good", "Fair", "Poor"] DAY_PROBS = [0.35, 0.45, 0.20] # if sum(DAY_PROBS) != 1: raise("Probabilities must sum to 1.") #...But computer numbers may not add to exactly 1: give it some slack: if abs(sum(DAY_PROBS)-1)>10^(-12): raise("Probabilities must sum to 1.") if min(DAY_PROBS) < 0: raise("Probabilities cannot be negative.") DEMAND_VALS = list(range(40,101,10)) DEMAND_PROBS = [[0.03,0.05,0.15,0.20,0.35,0.15,0.07], [0.10,0.18,0.40,0.20,0.08,0.04,0.00], [0.44,0.22,0.16,0.12,0.06,0.00,0.00]]

def day_type_rv(n_types = len(G.DAY_LIST), type_probs = G.DAY_PROBS): """ Random draw (with given probabilities) for the type of day. Returns an integer from 0 to (n_types-1) """ return rv_discrete(values=(range(n_types),type_probs)).rvs(size=1)[0]

def demand_rv(type_nbr, demand_vals = G.DEMAND_VALS, demand_probs = G.DEMAND_PROBS): """ Random draw (with given probabilities) for the type of day. Returns an integer from 0 to (n_types-1) """ return rv_discrete( values=(demand_vals,demand_probs[type_nbr])).rvs(size=1)[0]

def sim_month(days=20,Q=50,verbose=False): """ Monte Carlo simulation of 'days' (default is 20 days) of newspaper sales. Returns a list of the daily results: 1. Day number 2. Type of newsday (string) 3. Daily demand for newspapers 4. Daily revenue 5. Daily lost profit due to unsold papers 6. Daily scrap value 7. Daily profit """ dailies = [] #holds the daily results for i in range(days): day_type_nbr = day_type_rv() #random draw for the type of day demand = demand_rv(day_type_nbr) #random draw for demand based on the type of day revenue = min(demand,Q)*G.PRICE #revenue from direct sales cost = Q*G.COST #cost to buy Q papers lost_profit = max(0,demand-Q)*G.ILLWILL #penalty for unmet demand salvage = max(0,(Q-demand))*G.SALVAGE #salvage value for unsold papers daily_profit = revenue - cost - lost_profit + salvage #net profit dailies.append([i,G.DAY_LIST[day_type_nbr],demand,revenue, lost_profit,salvage,daily_profit]) if verbose: print("Day %d: %s, demand = %d, rev = %5.2f, lost = %5.2f, scrap = %5.2f, profit = %5.2f" % (i,G.DAY_LIST[day_type_nbr],demand, revenue,lost_profit,salvage,daily_profit)) return dailies

# returns confidence interval of mean def confIntMean(a, conf=0.95): """ Returns the lower and upper confidence interval limits for the data in list "a". Default is 95% CI. Assumes that the data comes from a normal distribution (that assumption is not so important if the number of observations is larger than 30 observations). """ mean, sem, m = np.mean(a), st.sem(a), st.t.ppf((1+conf)/2., len(a)-1) return mean - m*sem, mean + m*sem

np.random.seed(12345)

""" How to time a function: import timeit print(timeit.timeit("sim_month(Q=60)",number=100,setup="from __main__ import sim_month")) """

def sim_trials(reps=400,days=20,qty=60,conf_lvl=0.95, verbose=False,plotit=False,feedback=False): """ Run 'reps' number of sim_month trials with given number of 'days' in each month qty = order quantity of papers each day conf_lvl = confidence level for CIs If verbose is True, print interim counts If plotit is True, plot a histogram of the monthly profits for given Q Return average monthly profits for the given Q """ profits = [] # List to hold the monthly profit values #clock_it() # Start of the timing period for trial in range(reps): daily_vals = sim_month(days=days,Q=qty) monthly_profit = sum([row[6] for row in daily_vals]) # Sum daily profits profits.append(monthly_profit) if verbose & (trial%10 == 0): print("%d..." % trial,end="", flush=True) # Don't buffer: show it now!) if verbose & (trial%100 == 0): print() # Don't cram too much on one line... #print("It took %f seconds to run %d trials." % (clock_it(),REPS)) if verbose: print() print("For %d replications of %d days each, " % (reps,days) ) print("the average monthly profit was $%6.2f." % np.mean(profits)) cl_pct = conf_lvl*100. print("The %5.3f%% confidence interval " % cl_pct) print("for monthly profit is: (%6.2f, %6.2f)" % confIntMean(profits)) # Plot histogram of the total profits of all the trials plt.clf() # Clear the figure (i.e., plot) plt.hist(profits) # plot a histogram plt.title("Total Profits (%d Days, Q=%d)" % (days,qty)) # Add a title... plt.xlabel("Value") # and an x-axis label... plt.ylabel("Frequency") # and a y-axis label... plt.show() # Now, show the plot! return profits

print() print("News Dealer Problem") Q_range = range(40,101,10) profit_results = [] for Q in Q_range: print("For order quantity Q = %d:" % Q) profits = sim_trials(qty=Q) avg_profit = np.mean(profits) lower_CI, upper_CI = confIntMean(profits) profit_results.append((np.mean(profits),confIntMean(profits), (min(profits),max(profits)))) print("-----------------------------------------------------------") print()

avg_profits = [row[0] for row in profit_results] lower_CIs = [row[1][0] for row in profit_results] upper_CIs = [row[1][1] for row in profit_results] min_vals = [row[2][0] for row in profit_results] max_vals = [row[2][1] for row in profit_results]

print(" Q Mean CI Lower CI Upper CI Upper CI Upper") print("----- -------- -------- -------- -------- --------") i=0 for Q in Q_range: print("%3d %8.2f %8.2f %8.2f %8.2f %8.2f" % (Q,avg_profits[i],lower_CIs[i],upper_CIs[i], min_vals[i],max_vals[i])) i += 1

plt.clf() plt.plot(list(Q_range),avg_profits,"r-^") plt.vlines(list(Q_range),min_vals,max_vals,"r",linewidth=2) plt.title("Average Profits vs Daily Order Qty") plt.xlabel("Number of Papers Ordered") plt.ylabel("Average Profit")

What do I need to modify to fit the scenario?

1) the expected number of customers per day is 10.2; (2) the expected number of bagels is 2 dozen per customer; and (3) the expected number of bagels needed per day (multiply the answers from (1) and (2)) is 20.4. Based on those values, the owner would make 20 dozen bagels each day, with (4) expected cost of 20*$5.80 = $116 per day, (5) expected income of 20*$8.40 = 168, and (6) an expected profit of $168-$116 = $52