Jaką funkcję mogę użyć w Pythonie, jeśli chcę wypróbować skrócone prawo zasilania sieciowego?Próbka skróconego prawa zasilania liczbami całkowitymi w języku Python?
To znaczy, biorąc pod uwagę dwa parametry a i m wygenerować losową liczbę całkowitą w zakresie x[1,m) który następuje proporcjonalny dystrybucyjnego do 1/x^a.
Szukałem około numpy.random, ale nie znalazłem tej dystrybucji.
AFAIK, ani NumPy ani Scipy nie określają tego rozkładu dla ciebie. Jednak użycie scipy łatwo jest zdefiniować własną funkcję rozkładu dyskretnego użyciu 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()
Wygląda na to, że integrujesz pmf tak, jakby był ciągły, i biorąc obszar od 1 do 2, aby wymyślić p (1), między 2 a 3 dla p (2), itp., Czy to prawda? Jeśli tak, na przykład myślę, że musisz emulować Spinal Tap i przejść do 11, aby uzyskać p (10). Twoje 'const' będzie dostosowane poprzez dodanie' (m + 1) ** k' w mianowniku. Czy nie rozumiem? – pjs
@pjs: Biorę pdf jako * ciągłą * funkcję '1/x ** a'. Więc nie ma żadnej integracji w odstępach [1,2], [2,3] itd. Jednak zintegrowałem (ręcznie), aby znaleźć formuły dla 'const' i' _ppf', odwrotność 'cdf' . Myślę, że * dobrze to zrozumiałem, ale mogę się mylić. (Próbowałem twojej sugestii, ale przesuwa ona domenę na '[1, 11]', więc jeśli dobrze cię rozumiem, to nie przechodzi to podstawowej kontroli poprawności.) Przy okazji, czym jest Spinal Tap odnoszące się do tutaj? – unutbu
Spinal Tap był mockumentarnym filmem o heavy metalowym zespole. Wyróżnili się oni z innych zespołów, mając wzmacniacze do 11. – pjs
nie używam Python, więc zamiast błędy składniowe ryzyka postaram się opisać rozwiązanie algorytmicznie. Jest to inwersja dyskretna o charakterze brutalnej siły. Powinno to dość łatwo przetłumaczyć na język Python. Zakładam indeksowanie 0 dla tablicy.
Ustawienia:
generowanie tablicy cdf wielkości m z cdf[0] = 1 jako pierwszy wpis, cdf[i] = cdf[i-1] + 1/(i+1)**a dla pozostałych pozycji.
Skaluj wszystkie wpisy, dzieląc cdf[m-1] na każdy - teraz są to wartości CDF.
Zastosowanie:
cdf[] aż znajdziesz wpis większa niż jednolite . Zwróć indeks + 1 jako swoją wartość x.Powtórz czynność dla dowolnej liczby x-wartości.
Na przykład, z a,m = 2,10, obliczyć prawdopodobieństwo bezpośrednio jako:
[0.6452579827864142, 0.16131449569660355, 0.07169533142071269, 0.04032862392415089, 0.02581031931145657, 0.017923832855178172, 0.013168530260947229, 0.010082155981037722, 0.007966147935634743, 0.006452579827864143]
i CDF jest:
[0.6452579827864142, 0.8065724784830177, 0.8782678099037304, 0.9185964338278814, 0.944406753139338, 0.9623305859945162, 0.9754991162554634, 0.985581272236501, 0.9935474201721358, 1.0]
Podczas generowania, jeżeli mam Uniform wynik 0,90 wrócę x=4 ponieważ 0.918 ... to pierwszy wpis w pliku CDF większy niż mój jednolity.
Jeśli martwisz się szybkością, możesz zbudować tablicę aliasów, ale z rozkładem geometrycznym prawdopodobieństwo wcześniejszego zakończenia wyszukiwania liniowego w tablicy jest dość wysokie. W podanym przykładzie, na przykład, zakończysz pierwszy podgląd prawie 2/3 czasu.
Dlaczego po prostu nie próbować odrzucania z wbudowanymi dystrybucjami prawa mocy? –