Simulazione Montecarlo: una guida pratica

In questo esempio assumiamo di voler assemblare tre blocchi all’interno di un contenitore di una determinata larghezza. Le dimensioni nominali sono riportate nell’immagine sottostante. Vediamo che da progetto c’è uno spazio nominale di 0,5 mm.

Tuttavia, le dimensioni reali di questi tre blocchi e di un contenitore possono variare per motivi tecnici. A scopo dimostrativo, assumiamo che nessuna di queste variazioni segua la distribuzione normale. I tre blocchi seguiranno le distribuzioni triangolari illustrate di seguito e le dimensioni del contenitore seguiranno una distribuzione uniforme in un intervallo di +/-0,1 mm.

Ora, calcolando semplicemente i valori estremi, possiamo vedere che nello scenario peggiore i blocchi hanno 17 mm e un contenitore ha una larghezza di soli 16,4 mm, il che significa che, in questo caso, non possiamo inserirli tutti insieme.

La domanda è: qual è la probabilità di non riuscire a inserire tutti i blocchi in un contenitore?

In questo caso, le relazioni tra i blocchi appaiono come = Y = X4 – X1 – X2 – X3

Modificando il codice precedente otteniamo una funzione per campionare la distribuzione triangolare. Lo stesso possiamo farlo per ottenere la funzione di campionamento della distribuzione uniforme.

def mc_triangle(left, mode, right, samples):
results = []
for _ in range(samples):
results.append(rnd.triangular(left, mode, right))
return np.array(results)

Un codice di base modificato per la simulazione MC:

# configuration
s = 100000 # number of samples
limit = 0 # limit from specification
# components
component_1 = mc_triangle(3.8, 4.0, 4.4,s)
component_2 = mc_triangle(5.9, 6.0, 6.5,s)
component_3 = mc_triangle(5.8, 6.0, 6.1,s)
component_4 = mc_uniform(16.4, 16.6, s)
# relationships
total = component_4 - component_1 - component_2 - component_3
# success conditions
probability = np.sum(total < limit)/len(total)*100
print("Probability of not fitting the blocks: ", round(probability, 4), "%")

Dopo aver eseguito il codice precedente otteniamo una risposta nell’ordine di:

Probability of not fitting the blocks:  5.601 %

Dopo aver controllato la media e la deviazione standard, possiamo dire ancora di più sulla dimensione del divario atteso:

print('The mean gap width will be {}mm with a standard deviation of {}mm.'.format(round(np.mean(total), 2) ,round(np.std(total),2)))

The mean gap width will be 0.33mm with a standard deviation of 0.2mm.

Ora plottiamo l’istogramma del parametro stimato (larghezza del gap) per mostrare che segue una distribuzione normale nonostante nessuna delle distribuzioni in input era di questo tipo.

Ci si può ora chiedere come varierà il risultato all’aumentare del numero di campioni. Per verificarlo, si veda il grafico seguente che mostra la stima della larghezza del gap con un intervallo di confidenza del 95%, a seconda della dimensione del campione (da 100 a 7000 campioni):

Da questo grafico è evidente che la media del valore stimato non cambia in modo significativo, ma lo diffusione diminuisce con il numero di campioni. Ciò significa che con una nuova esecuzione della simulazione, campioni più grandi consentono di ottenere una minore diffusione dei risultati.

Tuttavia, a un certo punto l’aggiunta di più campioni non è più utile.