Generieren Sie Zufallszahlen aus der "geneigten Gleichverteilung" aus der mathematischen Theorie

9

Für einen bestimmten Zweck muss ich Zufallszahlen (Daten) aus einer "geneigten gleichmäßigen" Verteilung generieren. Die "Steigung" dieser Verteilung kann in einem angemessenen Intervall variieren, und dann sollte sich meine Verteilung basierend auf der Steigung von gleichmäßig zu dreieckig ändern. Hier ist meine Ableitung:

Geben Sie hier die Bildbeschreibung ein

Machen wir es einfach und generieren Daten von bis (blau, rot ist gleichmäßige Verteilung). Um die Wahrscheinlichkeitsdichtefunktion der blauen Linie zu erhalten, benötige ich nur die Gleichung dieser Linie. Somit:B.0B

f(x)=tg(φ)x+Y(0)

und da (Bild):

tg(φ)=1/BY(0)B/2Y(0)=1Btg(φ)B2

Wir haben das:

f(x)=tg(φ)x+(1Btg(φ)B2)

Da PDF ist, ist CDF gleich:f(x)

F(x)=tg(φ)x22+x(1Btg(φ)B2)

Jetzt machen wir einen Datengenerator. Die Idee ist, dass, wenn ich repariere , Zufallszahlen berechnet werden können, wenn ich Zahlen aus aus einer gleichmäßigen Verteilung wie hier beschrieben erhalte . Wenn ich also 100 Zufallszahl aus meiner Verteilung muß mit festen , dann für jede aus gleichmäßiger Verteilung gibt es von „geneigter Verteilung“, und kann wie folgt berechnet werden:xφ,Bxφ , B t i ((0,1)φ,Btix i x(0,1)xix

tg(φ)xi22+xi(1Btg(φ)B2)ti=0

Aus dieser Theorie habe ich Code in Python erstellt, der wie folgt aussieht:

import numpy as np
import math
import random
def tan_choice():
    x = random.uniform(-math.pi/3, math.pi/3)
    tan = math.tan(x)
    return tan

def rand_shape_unif(N, B, tg_fi):
    res = []
    n = 0
    while N > n:
        c = random.uniform(0,1)
        a = tg_fi/2
        b = 1/B - (tg_fi*B)/2
        quadratic = np.poly1d([a,b,-c])
        rots = quadratic.roots
        rot = rots[(rots.imag == 0) & (rots.real >= 0) & (rots.real <= B)].real
        rot = float(rot)
        res.append(rot)
        n += 1
    return res

def rand_numb(N_, B_):
    tan_ = tan_choice()
    res = rand_shape_unif(N_, B_, tan_)
    return res

Aber die Zahlen, aus rand_numbdenen generiert wird, sind sehr nahe bei Null oder bei B (was ich als 25 eingestellt habe). Es gibt keine Varianz, wenn ich 100 Zahlen generiere, sind alle nahe bei 25 oder alle nahe bei Null. In einem Lauf:

num = rand_numb(100, 25)
numb
Out[140]: 
[0.1063241766836174,
 0.011086243095907753,
 0.05690217839063588,
 0.08551031241199764,
 0.03411227661295121,
 0.10927087752739746,
 0.1173334720516189,
 0.14160616846114774,
 0.020124543145515768,
 0.10794924067959207]

In meinem Code muss also etwas sehr falsch sein. Kann mir jemand bei meiner Ableitung oder meinem Code helfen? Ich bin jetzt verrückt danach, ich kann keinen Fehler sehen. Ich nehme an, dass R-Code mir ähnliche Ergebnisse liefert.

Robert
quelle
2
Wenn Sie nur Zufallszahlen generieren müssen, müssen Sie die Verteilung überhaupt nicht berechnen. Wirf einfach Pfeile auf dein Bild und behalte ihre x-Koordinaten bei. Wenn ein Pfeil jedoch im linken Dreieck mit der Bezeichnung " " landet , ändere seine x-Koordinate von in . Zum Beispiel gibt alle Werte und (einen reeller Parameter, die bei gegebenen Werten zwischen und , Ihre Verteilungen erzeugen) und Satz die Anzahl von Zufallswerten sein , das Sie benötigen. Hier ist Code:x B - x - 1 1ϕxBxBtheta11nRx<-runif(n,-1,1);x<-(ifelse(runif(n,-1,1)>theta*x,-x,x)+1)*(B/2)
whuber

Antworten:

9

Ihre Ableitung ist in Ordnung. Beachten Sie, dass Sie einschränken müssen , um eine positive Dichte für zu erhalten In Ihrem Code sollten Sie also zwischen , hier schlägt Ihr Code fehl.B 2 tan ϕ < 2. B = 25 ϕ ± tan - 1 2(0,B)

B2tanϕ<2.
B=25ϕ±tan12625

Sie können (und sollten) die Verwendung eines quadratischen Lösers vermeiden und dann die Wurzeln zwischen 0 und auswählen . Die zu lösende quadratische Polynomgleichung in ist mit Durch Konstruktion ist und ; auch steigt auf .x F ( x ) = t F ( x ) = 1Bx

F(x)=t
F(x)=12tanϕx2+(1BB2tanϕ)x.
F(0)=0F(B)=1F(0,B)

Daraus ist leicht ersichtlich, dass, wenn , der Teil der Parabel, an dem wir interessiert sind, ein Teil der rechten Seite der Parabel ist und die zu behaltende Wurzel die höchste der beiden Wurzeln ist ist Im Gegenteil, wenn , steht die Parabel auf dem Kopf und wir sind an ihrer linken interessiert Teil. Die zu behaltende Wurzel ist die niedrigste. Unter Berücksichtigung des Vorzeichens von scheint es, dass dies dieselbe Wurzel ist (dh die mit ) wie im ersten Fall.tanϕ>0

x=1tanϕ(B2tanϕ1B+(B2tanϕ1B)2+2tanϕt.)
tanϕ<0tanϕ+Δ

Hier ist ein R-Code.

phi <- pi/8; B <- 2
f <- function(t) (-(1/B - 0.5*B*tan(phi)) + 
       sqrt( (1/B - 0.5*B*tan(phi))**2 + 2 * tan(phi) * t))/tan(phi)
hist(f(runif(1e6)))

Histogramm 1

Und mit :ϕ<0

phi <- -pi/8
hist(f(runif(1e6)))

Geben Sie hier die Bildbeschreibung ein

Elvis
quelle
1
Ich habe einen Fehler gemacht, weil ich meinen Winkel außerhalb der Grenzen gesetzt habe, ich verstehe. Aber Ihre Erklärung, warum ich die Verwendung des numerischen Lösers von vermeiden sollte, ist für mich immer noch neblig. Können Sie bitte versuchen, es näher zu erklären? Ich liebe es, es zu bekommen. F(x)
Robert
@ Robert Ich denke, dein Code funktioniert gut, wenn der Wert von korrekt ist. Es verhindert jedoch, dass Sie potenzielle Probleme erkennen (was ist, wenn keine Lösung zwischen 0 und ? Oder wenn beide Lösungen vorliegen? Oder wenn es keine echte Lösung gibt?). Die zusätzliche Arbeit, um die Verwendung des vorgefertigten Lösers zu vermeiden, lohnt sich. B.ϕB
Elvis