Quale funzione posso usare in Python se voglio campionare una legge di potenza a tronchi interi?Esempio di legge di potenza intera troncata in Python?
Cioè, dato due parametri a e m, generare un intero casuale x nell'intervallo [1,m) che segue una distribuzione proporzionale 1/x^a.
Ho cercato circa numpy.random, ma non ho trovato questa distribuzione.
AFAIK, né NumPy né Scipy definiscono questa distribuzione per voi. Tuttavia, utilizzando SciPy è facile definire la propria funzione distribuzione discreta con scipy.rv_discrete:
import numpy as np
import scipy.stats as stats
import matplotlib.pyplot as plt
def truncated_power_law(a, m):
x = np.arange(1, m+1, dtype='float')
pmf = 1/x**a
pmf /= pmf.sum()
return stats.rv_discrete(values=(range(1, m+1), pmf))
a, m = 2, 10
d = truncated_power_law(a=a, m=m)
N = 10**4
sample = d.rvs(size=N)
plt.hist(sample, bins=np.arange(m)+0.5)
plt.show()
Sembra che tu stia integrando il pmf come se fosse continuo, e prendendo l'area tra 1 e 2 per ottenere p (1), tra 2 e 3 per p (2), ecc., È giusto? Se è così, per il tuo esempio penso che devi emulare Spinal Tap e andare su 11 per ottenere p (10). Il tuo 'const' sarebbe aggiustato avendo' (m + 1) ** k' nel denominatore. O sto fraintendendo? – pjs
@pjs: sto prendendo il pdf per la * funzione * continua '1/x ** a'. Quindi non c'è integrazione su intervalli [1,2], [2,3], ecc. Tuttavia, ho integrato (a mano) per trovare le formule per 'const' e' _ppf', l'inverso del 'cdf' . Penso * Ho capito bene, ma potrei sbagliarmi. (Ho provato il tuo suggerimento, ma sposta il dominio a '[1, 11]', quindi se sto capendo correttamente, questo non passa un controllo di base sulla sanità mentale.) A proposito, che cos'è Spinal Tap che si riferisce a Qui? – unutbu
Spinal Tap era un film mockumentary su una band heavy metal. Si sono distinti dalle altre band avendo i loro amplificatori andare a 11. – pjs
Non faccio uso di Python, quindi piuttosto che errori di sintassi di rischio Cercherò di descrivere la soluzione algoritmicamente. Questa è un'inversione discreta a forza bruta. Dovrebbe tradurre abbastanza facilmente in Python. Sto assumendo l'indicizzazione basata su 0 per l'array.
Setup:
generare una serie di dimensioni cdfm con cdf[0] = 1 come prima voce, cdf[i] = cdf[i-1] + 1/(i+1)**a per le voci rimanenti.
Ridimensionare tutte le voci dividendo cdf[m-1] in ciascuna - ora sono effettivamente valori CDF.
Usage:
cdf[] fino a trovare una voce più grande del tuo uniforme. Restituisce l'indice + 1 come valore x.Ripetere per il numero di valori x desiderato.
Ad esempio, con a,m = 2,10, a calcolare le probabilità direttamente come:
[0.6452579827864142, 0.16131449569660355, 0.07169533142071269, 0.04032862392415089, 0.02581031931145657, 0.017923832855178172, 0.013168530260947229, 0.010082155981037722, 0.007966147935634743, 0.006452579827864143]
e CDF è:
[0.6452579827864142, 0.8065724784830177, 0.8782678099037304, 0.9185964338278814, 0.944406753139338, 0.9623305859945162, 0.9754991162554634, 0.985581272236501, 0.9935474201721358, 1.0]
Quando si generano, se ho ottenuto un risultato uniforme di 0,90 mi di ritorno x=4 perché 0.918 ... è la prima voce CDF più grande della mia uniforme.
Se si è preoccupati della velocità, è possibile creare una tabella alias, ma con un decadimento geometrico la probabilità di terminazione anticipata di una ricerca lineare attraverso l'array è piuttosto elevata. Con l'esempio dato, ad esempio, terminerai la prima volta circa i 2/3 del tempo.
Perché non fare solo il campionamento del rifiuto con le distribuzioni di legge di potenza incorporate? –