Simulazione Monte Carlo con Python: creazione di un istogramma al volo

Question

Ho una domanda concettuale sulla creazione di un istogramma al volo con Python. Sto cercando di capire se c'è un buon algoritmo o forse un pacchetto esistente.

Ho scritto una funzione, che esegue una simulazione Monte Carlo, viene chiamata 1.000.000.000 di volte e restituisce un numero mobile a 64 bit alla fine di ogni analisi. Sotto è detta funzione:

def MonteCarlo(df,head,span): 
    # Pick initial truck 
    rnd_truck = np.random.randint(0,len(df)) 
    full_length = df['length'][rnd_truck] 
    full_weight = df['gvw'][rnd_truck] 

    # Loop using other random trucks until the bridge is full 
    while True: 
     rnd_truck = np.random.randint(0,len(df)) 
     full_length += head + df['length'][rnd_truck] 
     if full_length > span: 
      break 
     else: 
      full_weight += df['gvw'][rnd_truck] 

    # Return average weight per feet on the bridge 
    return(full_weight/span) 

df è un oggetto dataframe Pandas avente colonne etichettate come 'length' e 'gvw', che sono lunghezze camion e pesi, rispettivamente. head è la distanza tra due camion consecutivi, span è la lunghezza del ponte. La funzione posiziona casualmente i camion sul ponte finché la lunghezza totale del convoglio è inferiore alla lunghezza del ponte. Infine, calcola il peso medio dei camion esistenti sul ponte per piede (peso totale esistente sul ponte diviso per la lunghezza del ponte).

Come risultato, vorrei creare un istogramma tabulare che mostri la distribuzione dei valori restituiti, che possono essere tracciati successivamente. Ho avuto alcune idee in mente:

1. Continuare a raccogliere i valori restituiti in un vettore NumPy, quindi utilizzare le funzioni di istogramma esistenti una volta che l'analisi Montecarlo è completata. Questo non sarebbe fattibile, poiché se il mio calcolo fosse corretto, avrei bisogno di 7,5 GB di memoria solo per quel vettore (1.000.000.000 64 bit float ~ 7.5 GB)

2. Inizializzare una matrice numpy con un intervallo e un numero di contenitori determinati . Aumentare il numero di elementi nel contenitore corrispondente di uno alla fine di ogni corsa. Il problema è che non conosco l'intervallo di valori che otterrei. L'impostazione di un istogramma con un intervallo e una dimensione bin appropriata non è nota. Devo anche capire come assegnare valori ai bin corretti, ma penso che sia fattibile.

3. Fallo in qualche modo al volo. Modificare intervalli e dimensioni dei contenitori ogni volta che la funzione restituisce un numero. Sarebbe troppo complicato scrivere da zero, credo.

Beh, scommetto che potrebbe esserci un modo migliore per gestire questo problema. Qualche idea sarebbe gradita!

In una seconda nota, ho provato a eseguire la funzione sopra per 1.000.000.000 di volte solo per ottenere il valore più grande che è calcolato (il frammento di codice è sotto). E questo richiede circa un'ora quando span = 200. Il tempo di calcolo aumenterebbe se lo eseguivo per intervalli più lunghi (il ciclo while è più lungo per riempire il ponte con i camion). C'è un modo per ottimizzare ciò che pensi?

max_w = 0 
i = 1 
    while i < 1000000000: 
     if max_w < MonteCarlo(df_basic, 15., 200.): 
      max_w = MonteCarlo(df_basic, 15., 200.) 
    i += 1 
print max_w 

Grazie!

Answer 1

Ecco una possibile soluzione, con dimensione del contenitore fissa e bin della forma [k * dimensione, (k + 1) * dimensione [. La funzione finalizebins restituisce due liste: una con bin count (a), e l'altra (b) con bin lower bounds (il limite superiore viene dedotto aggiungendo binsize).

import math, random 

def updatebins(bins, binsize, x): 
    i = math.floor(x/binsize) 
    if i in bins: 
     bins[i] += 1 
    else: 
     bins[i] = 1 

def finalizebins(bins, binsize): 
    imin = min(bins.keys()) 
    imax = max(bins.keys()) 
    a = [0] * (imax - imin + 1) 
    b = [binsize * k for k in range(imin, imax + 1)] 
    for i in range(imin, imax + 1): 
     if i in bins: 
      a[i - imin] = bins[i] 
    return a, b 

# A test with a mixture of gaussian distributions 

def check(n): 
    bins = {} 
    binsize = 5.0 
    for i in range(n): 
     if random.random() > 0.5: 
      x = random.gauss(100, 50) 
     else: 
      x = random.gauss(-200, 150) 
     updatebins(bins, binsize, x) 
    return finalizebins(bins, binsize) 

a, b = check(10000) 

# This must be 10000 
sum(a) 

# Plot the data 
from matplotlib.pyplot import * 
bar(b,a) 
show()