Generuj liczby losowe o podanym (liczbowym) rozkładzie

Mam plik z pewnymi prawdopodobieństwami dla różnych wartości np:

1 0.1
2 0.05
3 0.05
4 0.2
5 0.4
6 0.2

Chciałbym wygenerować liczby losowe za pomocą tej dystrybucji. Czy istnieje moduł, który to obsługuje? Samodzielne kodowanie jest dość proste (zbuduj funkcję gęstości kumulacyjnej, wygeneruj losową wartość [0,1] i wybierz odpowiednią wartość), ale wydaje się, że to powinien być powszechny problem i prawdopodobnie ktoś stworzył funkcję / moduł to.

Potrzebuję tego, ponieważ chcę wygenerować listę urodzin (które nie są zgodne z żadną dystrybucją w randommodule standardowym ).

scipy.stats.rv_discretemoże być tym, czego chcesz. Możesz podać swoje prawdopodobieństwa za pomocą valuesparametru. Następnie możesz użyć rvs()metody obiektu dystrybucji, aby wygenerować liczby losowe.

Jak zauważył w komentarzach Eugene Pakhomov, można również przekazać pparametr słowa kluczowego numpy.random.choice()np

numpy.random.choice(numpy.arange(1, 7), p=[0.1, 0.05, 0.05, 0.2, 0.4, 0.2])

Jeśli korzystasz z Pythona 3.6 lub nowszego, możesz korzystać random.choices()z biblioteki standardowej - zobacz odpowiedź Marka Dickinsona .

Począwszy od Pythona 3.6, istnieje rozwiązanie tego problemu w standardowej bibliotece Pythona, a mianowicie random.choices.

Przykładowe użycie: skonfigurujmy populację i wagi odpowiadające tym w pytaniu PO:

>>> from random import choices
>>> population = [1, 2, 3, 4, 5, 6]
>>> weights = [0.1, 0.05, 0.05, 0.2, 0.4, 0.2]

Teraz choices(population, weights)generuje pojedynczą próbkę:

>>> choices(population, weights)
4

Opcjonalny argument zawierający tylko słowo kluczowe kpozwala zażądać więcej niż jednej próbki na raz. Jest to cenne, ponieważ random.choicesprzed wygenerowaniem jakichkolwiek próbek należy wykonać pewne prace przygotowawcze przy każdym wywołaniu; generując wiele próbek na raz, musimy wykonać tę pracę przygotowawczą tylko raz. Tutaj generujemy milion próbek i używamy collections.Counterdo sprawdzenia, czy otrzymany rozkład z grubsza odpowiada podanym przez nas wagom.

>>> million_samples = choices(population, weights, k=10**6)
>>> from collections import Counter
>>> Counter(million_samples)
Counter({5: 399616, 6: 200387, 4: 200117, 1: 99636, 3: 50219, 2: 50025})

Zaletą generowania listy przy użyciu CDF jest możliwość korzystania z wyszukiwania binarnego. Chociaż potrzebujesz O (n) czasu i miejsca na przetwarzanie wstępne, możesz uzyskać k liczb w O (k log n). Ponieważ zwykłe listy Pythona są nieefektywne, możesz użyć arraymodule.

Jeśli nalegasz na stałą przestrzeń, możesz wykonać następujące czynności; O (n) czas, O (1) przestrzeń.

def random_distr(l):
    r = random.uniform(0, 1)
    s = 0
    for item, prob in l:
        s += prob
        if s >= r:
            return item
    return item  # Might occur because of floating point inaccuracies

Może jest już trochę późno. Ale możesz użyć numpy.random.choice(), przekazując pparametr:

val = numpy.random.choice(numpy.arange(1, 7), p=[0.1, 0.05, 0.05, 0.2, 0.4, 0.2])

(OK, wiem, że prosisz o folię termokurczliwą, ale może te domowe rozwiązania nie były wystarczająco zwięzłe według twoich upodobań. :-)

pdf = [(1, 0.1), (2, 0.05), (3, 0.05), (4, 0.2), (5, 0.4), (6, 0.2)]
cdf = [(i, sum(p for j,p in pdf if j < i)) for i,_ in pdf]
R = max(i for r in [random.random()] for i,c in cdf if c <= r)

Pseudo-potwierdziłem, że to działa, patrząc na wynik tego wyrażenia:

sorted(max(i for r in [random.random()] for i,c in cdf if c <= r)
       for _ in range(1000))

Napisałem rozwiązanie do pobierania losowych próbek z niestandardowej ciągłej dystrybucji .

Potrzebowałem tego do podobnego przypadku użycia do twojego (tj. Generowania losowych dat z podanym rozkładem prawdopodobieństwa).

Potrzebujesz tylko funkcji random_custDisti linki samples=random_custDist(x0,x1,custDist=custDist,size=1000). Reszta to dekoracja ^^.

import numpy as np

#funtion
def random_custDist(x0,x1,custDist,size=None, nControl=10**6):
    #genearte a list of size random samples, obeying the distribution custDist
    #suggests random samples between x0 and x1 and accepts the suggestion with probability custDist(x)
    #custDist noes not need to be normalized. Add this condition to increase performance. 
    #Best performance for max_{x in [x0,x1]} custDist(x) = 1
    samples=[]
    nLoop=0
    while len(samples)<size and nLoop<nControl:
        x=np.random.uniform(low=x0,high=x1)
        prop=custDist(x)
        assert prop>=0 and prop<=1
        if np.random.uniform(low=0,high=1) <=prop:
            samples += [x]
        nLoop+=1
    return samples

#call
x0=2007
x1=2019
def custDist(x):
    if x<2010:
        return .3
    else:
        return (np.exp(x-2008)-1)/(np.exp(2019-2007)-1)
samples=random_custDist(x0,x1,custDist=custDist,size=1000)
print(samples)

#plot
import matplotlib.pyplot as plt
#hist
bins=np.linspace(x0,x1,int(x1-x0+1))
hist=np.histogram(samples, bins )[0]
hist=hist/np.sum(hist)
plt.bar( (bins[:-1]+bins[1:])/2, hist, width=.96, label='sample distribution')
#dist
grid=np.linspace(x0,x1,100)
discCustDist=np.array([custDist(x) for x in grid]) #distrete version
discCustDist*=1/(grid[1]-grid[0])/np.sum(discCustDist)
plt.plot(grid,discCustDist,label='custom distribustion (custDist)', color='C1', linewidth=4)
#decoration
plt.legend(loc=3,bbox_to_anchor=(1,0))
plt.show()

Wydajność tego rozwiązania jest na pewno możliwa do poprawy, ale ja wolę czytelność.

Zrób listę elementów na podstawie ich weights:

items = [1, 2, 3, 4, 5, 6]
probabilities= [0.1, 0.05, 0.05, 0.2, 0.4, 0.2]
# if the list of probs is normalized (sum(probs) == 1), omit this part
prob = sum(probabilities) # find sum of probs, to normalize them
c = (1.0)/prob # a multiplier to make a list of normalized probs
probabilities = map(lambda x: c*x, probabilities)
print probabilities

ml = max(probabilities, key=lambda x: len(str(x)) - str(x).find('.'))
ml = len(str(ml)) - str(ml).find('.') -1
amounts = [ int(x*(10**ml)) for x in probabilities]
itemsList = list()
for i in range(0, len(items)): # iterate through original items
  itemsList += items[i:i+1]*amounts[i]

# choose from itemsList randomly
print itemsList

Optymalizacją może być normalizacja kwot za pomocą największego wspólnego dzielnika, tak aby lista docelowa była mniejsza.

Także, to może być ciekawe.

Inna odpowiedź, chyba szybsza :)

distribution = [(1, 0.2), (2, 0.3), (3, 0.5)]  
# init distribution  
dlist = []  
sumchance = 0  
for value, chance in distribution:  
    sumchance += chance  
    dlist.append((value, sumchance))  
assert sumchance == 1.0 # not good assert because of float equality  

# get random value  
r = random.random()  
# for small distributions use lineair search  
if len(distribution) < 64: # don't know exact speed limit  
    for value, sumchance in dlist:  
        if r < sumchance:  
            return value  
else:  
    # else (not implemented) binary search algorithm  

from __future__ import division
import random
from collections import Counter


def num_gen(num_probs):
    # calculate minimum probability to normalize
    min_prob = min(prob for num, prob in num_probs)
    lst = []
    for num, prob in num_probs:
        # keep appending num to lst, proportional to its probability in the distribution
        for _ in range(int(prob/min_prob)):
            lst.append(num)
    # all elems in lst occur proportional to their distribution probablities
    while True:
        # pick a random index from lst
        ind = random.randint(0, len(lst)-1)
        yield lst[ind]

Weryfikacja:

gen = num_gen([(1, 0.1),
               (2, 0.05),
               (3, 0.05),
               (4, 0.2),
               (5, 0.4),
               (6, 0.2)])
lst = []
times = 10000
for _ in range(times):
    lst.append(next(gen))
# Verify the created distribution:
for item, count in Counter(lst).iteritems():
    print '%d has %f probability' % (item, count/times)

1 has 0.099737 probability
2 has 0.050022 probability
3 has 0.049996 probability 
4 has 0.200154 probability
5 has 0.399791 probability
6 has 0.200300 probability

bazując na innych rozwiązaniach, generujesz dystrybucję akumulacyjną (jako liczbę całkowitą lub zmiennoprzecinkową, jak chcesz), a następnie możesz użyć bisect, aby przyspieszyć

to jest prosty przykład (użyłem tutaj liczb całkowitych)

l=[(20, 'foo'), (60, 'banana'), (10, 'monkey'), (10, 'monkey2')]
def get_cdf(l):
    ret=[]
    c=0
    for i in l: c+=i[0]; ret.append((c, i[1]))
    return ret

def get_random_item(cdf):
    return cdf[bisect.bisect_left(cdf, (random.randint(0, cdf[-1][0]),))][1]

cdf=get_cdf(l)
for i in range(100): print get_random_item(cdf),

get_cdffunkcja będzie przekształcić od 20, 60, 10, 10 do 20, 20 + 60, 20 + 60 + 10 20 + 60 + 10 + 10

teraz wybieramy losową liczbę do 20 + 60 + 10 + 10 za pomocą, random.randinta następnie używamy połowy, aby szybko uzyskać rzeczywistą wartość

warto rzucić okiem na rozkłady próbkowania NumPy Random

Żadna z tych odpowiedzi nie jest szczególnie jasna ani prosta.

Oto jasna, prosta metoda, która na pewno zadziała.

umulate_normalize_probabilities pobiera słownik, pktóry odwzorowuje symbole na prawdopodobieństwa LUB częstotliwości. Wyświetla użyteczną listę krotek, z których można dokonać wyboru.

def accumulate_normalize_values(p):
        pi = p.items() if isinstance(p,dict) else p
        accum_pi = []
        accum = 0
        for i in pi:
                accum_pi.append((i[0],i[1]+accum))
                accum += i[1]
        if accum == 0:
                raise Exception( "You are about to explode the universe. Continue ? Y/N " )
        normed_a = []
        for a in accum_pi:
                normed_a.append((a[0],a[1]*1.0/accum))
        return normed_a

Plony:

>>> accumulate_normalize_values( { 'a': 100, 'b' : 300, 'c' : 400, 'd' : 200  } )
[('a', 0.1), ('c', 0.5), ('b', 0.8), ('d', 1.0)]

Dlaczego to działa

Etap akumulacji zamienia każdy symbol w przedział między nim a prawdopodobieństwem lub częstotliwością poprzednich symboli (lub 0 w przypadku pierwszego symbolu). Przedziały te mogą być używane do wybierania z (a tym samym próbkowania dostarczonego rozkładu), po prostu przechodząc przez listę, aż liczba losowa w przedziale 0,0 -> 1,0 (przygotowana wcześniej) będzie mniejsza lub równa punktowi końcowemu interwału bieżącego symbolu.

Normalizacja uwalnia nas od konieczności upewnić, że wszystko sum do pewnej wartości. Po normalizacji „wektor” prawdopodobieństw sumuje się do 1,0.

Reszta kodu dla selekcji i generowanie dowolnie długi próbki z rozkładu jest poniżej:

def select(symbol_intervals,random):
        print symbol_intervals,random
        i = 0
        while random > symbol_intervals[i][1]:
                i += 1
                if i >= len(symbol_intervals):
                        raise Exception( "What did you DO to that poor list?" )
        return symbol_intervals[i][0]


def gen_random(alphabet,length,probabilities=None):
        from random import random
        from itertools import repeat
        if probabilities is None:
                probabilities = dict(zip(alphabet,repeat(1.0)))
        elif len(probabilities) > 0 and isinstance(probabilities[0],(int,long,float)):
                probabilities = dict(zip(alphabet,probabilities)) #ordered
        usable_probabilities = accumulate_normalize_values(probabilities)
        gen = []
        while len(gen) < length:
                gen.append(select(usable_probabilities,random()))
        return gen

Stosowanie :

>>> gen_random (['a','b','c','d'],10,[100,300,400,200])
['d', 'b', 'b', 'a', 'c', 'c', 'b', 'c', 'c', 'c']   #<--- some of the time

Oto skuteczniejszy sposób na zrobienie tego:

Po prostu wywołaj następującą funkcję z tablicą „weights” (zakładając, że indeksy są odpowiednimi elementami) i nie. potrzebnych próbek. Funkcję tę można łatwo zmodyfikować w celu obsługi uporządkowanej pary.

Zwraca indeksy (lub elementy) próbkowane / pobierane (z wymianą) przy użyciu odpowiednich prawdopodobieństw:

def resample(weights, n):
    beta = 0

    # Caveat: Assign max weight to max*2 for best results
    max_w = max(weights)*2

    # Pick an item uniformly at random, to start with
    current_item = random.randint(0,n-1)
    result = []

    for i in range(n):
        beta += random.uniform(0,max_w)

        while weights[current_item] < beta:
            beta -= weights[current_item]
            current_item = (current_item + 1) % n   # cyclic
        else:
            result.append(current_item)
    return result

Krótka uwaga na temat koncepcji używanej w pętli while. Zmniejszamy wagę bieżącego przedmiotu ze skumulowanej beta, która jest skumulowaną wartością konstruowaną równomiernie losowo, i zwiększamy bieżący indeks w celu znalezienia przedmiotu, którego waga odpowiada wartości beta.