Generuj losowe liczby z „nachylonego rozkładu równomiernego” z teorii matematycznej

9

W jakimś celu muszę wygenerować losowe liczby (dane) z rozkładu „nachylonego równomiernie”. „Nachylenie” tego rozkładu może się zmieniać w rozsądnych odstępach czasu, a następnie mój rozkład powinien zmienić się z jednorodnego na trójkątny w zależności od nachylenia. Oto moje pochodzenie:

wprowadź opis zdjęcia tutaj

Uprośćmy to i wygeneruj dane od do (niebieski, czerwony to rozkład równomierny). Aby uzyskać funkcję gęstości prawdopodobieństwa niebieskiej linii, potrzebuję tylko równania tej linii. A zatem:0B

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

i ponieważ (zdjęcie):

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

Mamy to:

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

Ponieważ to PDF, CDF wynosi:f(x)

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

Teraz stwórzmy generator danych. Chodzi o to, że jeśli naprawię , liczby losowe można obliczyć, jeśli otrzymam liczby od z jednolitego rozkładu, jak opisano tutaj . Zatem jeśli potrzebuję 100 liczb losowych z mojego rozkładu ze stałym , to dla dowolnego z rozkładu jednolitego istnieje z „rozkładu nachylonego”, a można obliczyć jako:φ,Bx(0,1)φ,Bti(0,1)xix

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

Z tej teorii stworzyłem kod w Pythonie, który wygląda następująco:

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

Ale liczby generowane z rand_numbsą bardzo bliskie zeru lub B (które ustawiłem na 25). Nie ma wariancji, gdy generuję 100 liczb, wszystkie są bliskie 25 lub wszystkie są bliskie zeru. W jednym biegu:

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]

Więc musi być coś bardzo złego w moim kodzie. Czy ktoś może mi pomóc z moim wyprowadzeniem lub kodem? Jestem szalony z tego powodu, nie widzę żadnego błędu. Przypuszczam, że kod R da mi podobne wyniki.

Robert
źródło
2
Jeśli potrzebujesz wygenerować tylko liczby losowe, nie musisz w ogóle ustalać rozkładu. Po prostu rzuć rzutkami w swoje zdjęcie i zachowaj ich współrzędne x, ale kiedy strzałka wyląduje w lewym trójkącie oznaczonym „ ”, zmień współrzędną z na . Na przykład podaj dowolne wartości do i (prawdziwy parametr, który, gdy podano wartości od do , tworzy twoje rozkłady) i ustaw na liczbę potrzebnych losowych wartości. Oto kod:ϕxBxBtheta11nRx<-runif(n,-1,1);x<-(ifelse(runif(n,-1,1)>theta*x,-x,x)+1)*(B/2)
whuber

Odpowiedzi:

9

Twoje pochodzenie jest w porządku. Zauważ, że aby uzyskać dodatnią gęstość na , musisz ograniczyć W kodzie więc powinieneś wziąć między , tam właśnie zawiedzie Twój kod.(0,B)

B2tanϕ<2.
B=25ϕ±tan12625

Można (i powinno) Należy unikać używania kwadratowego Solver, a następnie wybierz korzenie między 0 i . Kwadratyczne równanie wielomianowe w do rozwiązania to z Z założenia i ; również wzrasta na .Bx

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

Z tego łatwo zauważyć, że jeśli , część paraboli, którą jesteśmy zainteresowani, jest częścią prawej strony paraboli, a korzeń do zachowania jest najwyższym z dwóch korzeni, że to Przeciwnie, jeśli , parabola jest do góry nogami, a my interesujemy się jej lewą stroną część. Korzeń do zachowania jest najniższy. Biorąc pod uwagę znak wydaje się, że jest to ten sam root (tj. Ten z ) niż w pierwszym przypadku.tanϕ>0

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

Oto trochę kodu R.

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

histogram 1

I z :ϕ<0

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

wprowadź opis zdjęcia tutaj

Elvis
źródło
1
Popełniłem błąd, ponieważ ustawiłem swój kąt poza granicami, rozumiem. Ale twoje wyjaśnienie, dlaczego powinienem unikać omawiania liczb za pomocą jest dla mnie nadal mgliste. Czy możesz spróbować wyjaśnić więcej? uwielbiam to robić. F(x)
Robert
@Robert Myślę, że twój kod działa dobrze, jeśli wartość jest poprawna. Jednak zapobiega uchwyceniu potencjalnych problemów (co, jeśli żadne rozwiązanie nie jest pomiędzy 0 a ? Lub jeśli oba rozwiązania są? Lub jeśli nie ma prawdziwego rozwiązania?). Warto dodatkowej pracy, aby uniknąć korzystania z gotowego solvera. ϕB
Elvis