Generieren Sie Zufallszahlen mit einer bestimmten (numerischen) Verteilung

Ich habe eine Datei mit einigen Wahrscheinlichkeiten für verschiedene Werte, z.

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

Ich möchte mit dieser Verteilung Zufallszahlen generieren. Gibt es ein vorhandenes Modul, das dies behandelt? Es ist ziemlich einfach, selbst zu codieren (die kumulative Dichtefunktion erstellen, einen zufälligen Wert [0,1] generieren und den entsprechenden Wert auswählen), aber es scheint, dass dies ein häufiges Problem sein sollte und wahrscheinlich jemand eine Funktion / ein Modul für erstellt hat es.

Ich brauche das, weil ich eine Liste von Geburtstagen erstellen möchte (die keiner Verteilung im Standardmodul folgen random).

scipy.stats.rv_discretekönnte sein, was Sie wollen. Sie können Ihre Wahrscheinlichkeiten über den valuesParameter angeben. Sie können dann die rvs()Methode des Verteilungsobjekts verwenden, um Zufallszahlen zu generieren.

Wie Eugene Pakhomov in den Kommentaren hervorhob, können Sie auch einen pSchlüsselwortparameter an numpy.random.choice()z

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

Wenn Sie Python 3.6 oder höher verwenden, können Sie es random.choices()aus der Standardbibliothek verwenden - siehe die Antwort von Mark Dickinson .

Seit Python 3.6 gibt es dafür eine Lösung in der Standardbibliothek von Python, nämlich random.choices.

Anwendungsbeispiel: Lassen Sie uns eine Grundgesamtheit und Gewichte einrichten, die denen in der Frage des OP entsprechen:

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

choices(population, weights)Generiert jetzt eine einzelne Stichprobe:

>>> choices(population, weights)
4

Mit dem optionalen Argument "Nur Schlüsselwörter" kkönnen mehrere Beispiele gleichzeitig angefordert werden. Dies ist wertvoll, da random.choicesjedes Mal, wenn es aufgerufen wird, einige Vorbereitungsarbeiten durchgeführt werden müssen, bevor Proben generiert werden. Durch die gleichzeitige Erzeugung vieler Proben müssen wir diese Vorbereitungsarbeiten nur einmal durchführen. Hier generieren wir eine Million Proben und collections.Counterüberprüfen, ob die Verteilung, die wir erhalten, ungefähr mit den von uns angegebenen Gewichten übereinstimmt.

>>> 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})

Ein Vorteil beim Generieren der Liste mit CDF besteht darin, dass Sie die binäre Suche verwenden können. Während Sie O (n) Zeit und Raum für die Vorverarbeitung benötigen, können Sie k Zahlen in O (k log n) erhalten. Da normale Python-Listen ineffizient sind, können Sie das arrayModul verwenden.

Wenn Sie auf konstantem Platz bestehen, können Sie Folgendes tun: O (n) Zeit, O (1) Raum.

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

Vielleicht ist es schon spät. Sie können jedoch Folgendes verwenden numpy.random.choice(), indem Sie den pParameter übergeben:

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

(OK, ich weiß, dass Sie nach Schrumpffolie fragen, aber vielleicht waren diese selbst entwickelten Lösungen einfach nicht prägnant genug für Ihren Geschmack. :-)

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)

Ich habe pseudo-bestätigt, dass dies funktioniert, indem ich die Ausgabe dieses Ausdrucks musterte:

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

Ich habe eine Lösung zum Zeichnen von Zufallsstichproben aus einer benutzerdefinierten kontinuierlichen Verteilung geschrieben .

Ich brauchte dies für einen ähnlichen Anwendungsfall wie Ihren (dh das Generieren von zufälligen Daten mit einer bestimmten Wahrscheinlichkeitsverteilung).

Sie brauchen nur die Funktion random_custDistund die Leitung samples=random_custDist(x0,x1,custDist=custDist,size=1000). Der Rest ist Dekoration ^^.

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()

Die Leistung dieser Lösung ist sicher verbesserungsfähig, aber ich bevorzuge die Lesbarkeit.

Erstellen Sie eine Liste mit Elementen, basierend auf 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

Eine Optimierung kann darin bestehen, Beträge durch den größten gemeinsamen Teiler zu normalisieren, um die Zielliste kleiner zu machen.

Auch das könnte interessant sein.

Eine andere Antwort, wahrscheinlich schneller :)

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]

Überprüfung:

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

Basierend auf anderen Lösungen generieren Sie eine kumulative Verteilung (als Ganzzahl oder Float, was auch immer Sie möchten). Anschließend können Sie die Halbierung verwenden, um sie schnell zu machen

Dies ist ein einfaches Beispiel (ich habe hier ganze Zahlen verwendet)

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),

Die get_cdfFunktion würde es von 20, 60, 10, 10 in 20, 20 + 60, 20 + 60 + 10, 20 + 60 + 10 + 10 konvertieren

Jetzt wählen wir eine Zufallszahl bis zu 20 + 60 + 10 + 10 aus und verwenden random.randintdann die Halbierung, um den tatsächlichen Wert schnell zu erhalten

Vielleicht möchten Sie einen Blick auf NumPy Random Sampling Distributionen werfen

Keine dieser Antworten ist besonders klar oder einfach.

Hier ist eine klare, einfache Methode, die garantiert funktioniert.

accumulate_normalize_probabilities verwendet ein Wörterbuch p, das Symbole Wahrscheinlichkeiten ODER Frequenzen zuordnet . Es gibt eine verwendbare Liste von Tupeln aus, aus denen ausgewählt werden kann.

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

Ausbeuten:

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

Warum es funktioniert

Die Akkumulation verwandelt jedes Symbol in ein Intervall zwischen sich und der Wahrscheinlichkeit oder Häufigkeit der vorherigen Symbole (oder 0 im Fall des ersten Symbols). Diese Intervalle können verwendet werden, um aus der Liste auszuwählen (und damit die bereitgestellte Verteilung abzutasten), indem Sie einfach durch die Liste gehen, bis die Zufallszahl in Intervall 0.0 -> 1.0 (zuvor vorbereitet) kleiner oder gleich dem Intervallendpunkt des aktuellen Symbols ist.

Das Normalisierung befreit uns von der Notwendigkeit, sicherzustellen, dass alles einen gewissen Wert hat. Nach der Normalisierung summiert sich der "Vektor" der Wahrscheinlichkeiten auf 1,0.

Der Rest des Codes zum Auswählen und Generieren einer beliebig langen Stichprobe aus der Verteilung ist unten aufgeführt:

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

Verwendung :

>>> 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

Hier ist ein effektiverer Weg , dies zu tun:

Rufen Sie einfach die folgende Funktion mit Ihrem 'Gewichte'-Array (unter der Annahme, dass die Indizes die entsprechenden Elemente sind) und der Nr. Auf. von Proben benötigt. Diese Funktion kann leicht geändert werden, um geordnete Paare zu handhaben.

Gibt Indizes (oder Artikel) zurück, die mit ihren jeweiligen Wahrscheinlichkeiten abgetastet / ausgewählt wurden (mit Ersatz):

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

Ein kurzer Hinweis zum Konzept der while-Schleife. Wir reduzieren das Gewicht des aktuellen Artikels aus dem kumulativen Beta, einem kumulativen Wert, der gleichmäßig zufällig erstellt wird, und erhöhen den aktuellen Index, um den Artikel zu finden, dessen Gewicht mit dem Beta-Wert übereinstimmt.