en_gros.py

import numpy as np
import matplotlib.pyplot as plt


def monte_carlo(n_trajects, n_clients, n_iter, capacity=None, plot=False):
	# maximum capacity by car is here unlimited
	if capacity is None:
		capacity = n_clients

	# draw within an uniform distribution
	draw = np.random.randint(0, n_trajects, (n_iter, n_clients))

	# count the number of car required
	n_car = np.zeros(n_iter)
	for i in range(n_iter):
		_, count = np.unique(draw[i], return_counts=True)
		n_car[i] = ((count - 1) // capacity + 1).sum()

	# compute the histogram
	unique, count = np.unique(n_car, return_counts=True)
	count = count / float(n_iter)

	# statistics
	mean = unique.mean()
	std = unique.std()

	if plot:
		plt.figure()
		plt.plot(unique, count, '.-')
		plt.title('One simulation (n_trajects = %d, n_clients = %d)'
				  % (n_trajects, n_clients))
		plt.xlabel('number of unique traject')
		plt.ylabel('probability distribution')

	return mean, std

n_trajects = 100  # number of trajects possible
n_iter = 100  # number of iteration for the simulation
n_stats = n_trajects * 3  # maximum number of clients
capacities = [1, 2, 3, 8]  # capacity of a car

# one simulation
monte_carlo(n_trajects, n_trajects // 2, n_iter, plot=False)

# multiple simulations for n_clients in range(n_stats)
plt.figure()
for capacity in capacities:
	means = np.zeros(n_stats)
	stds = np.zeros(n_stats)
	for i, n_clients in enumerate(range(n_stats)):
		means[i], stds[i] = monte_carlo(n_trajects, n_clients, n_iter,
						capacity=capacity, plot=False)

	plt.plot(range(n_stats), means, label=('capacity=%d' % capacity))
	plt.fill_between(range(n_stats), means - stds, means + stds,
				     facecolor='yellow', alpha=0.4)
	plt.title('Multiple simulations (n_trajects = %d)' % n_trajects)
	plt.xlabel('number of clients')
	plt.ylabel('number of car needed')
plt.legend(loc=0)
plt.grid()
plt.show()