Genere números aleatorios a partir de "distribución uniforme inclinada" a partir de la teoría matemática

9

Para algún propósito, necesito generar números aleatorios (datos) de la distribución "uniforme inclinado". La "pendiente" de esta distribución puede variar en un intervalo razonable, y luego mi distribución debería cambiar de uniforme a triangular en función de la pendiente. Aquí está mi derivación:

ingrese la descripción de la imagen aquí

Hagámoslo simple y generemos datos de a (azul, rojo es distribución uniforme). Para obtener la función de densidad de probabilidad de la línea azul solo necesito la ecuación de esa línea. Así:B0B

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

y desde (foto):

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

Tenemos que:

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

Como es PDF, CDF es igual a:f(x)

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

Ahora hagamos un generador de datos. La idea es que si soluciono , se pueden calcular números aleatorios si obtengo números de de una distribución uniforme como se describe aquí . Por lo tanto, si necesito 100 números aleatorios de mi distribución con fijo , entonces para cualquier de distribución uniforme hay de "distribución inclinada", y puede calcularse como:x ( 0 , 1 ) φ , B t i ( 0 , 1 )φ,Bx(0,1)φ,Bti(0,1) xxix

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

A partir de esta teoría, hice un código en Python que se parece a:

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

Pero los números generados rand_numbson muy cercanos a cero o B (que configuré como 25). No hay variación, cuando genero 100 números, todos ellos están cerca de 25 o todos están cerca de cero. En una carrera:

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]

Entonces debe haber algo muy mal en mi código. ¿Alguien puede ayudarme con mi derivación o código? Estoy loco por esto ahora, no puedo ver ningún error. Supongo que el código R me dará resultados similares.

Robert
fuente
2
Si solo necesita generar números aleatorios, no tiene que calcular la distribución en absoluto. Simplemente arroje dardos a su imagen y conserve sus coordenadas x, pero cuando un dardo aterrice en el triángulo izquierdo con la etiqueta " ", cambie su coordenada de a . Por ejemplo, proporcione cualquier valor a y (un parámetro real que, cuando se le dan valores entre y , produce sus distribuciones) y establezca el número de valores aleatorios que necesita. Aquí está el código:x B - x - 1 1ϕxBxBtheta11nRx<-runif(n,-1,1);x<-(ifelse(runif(n,-1,1)>theta*x,-x,x)+1)*(B/2)
whuber

Respuestas:

9

Tu derivación está bien. Tenga en cuenta que para obtener una densidad positiva en , debe restringir En su código , debe tomar entre , ahí es donde falla su código.(0,B)

B2tanϕ<2.
B=25ϕ±tan12625

Usted puede (y debe) evitar el uso de un programa de solución cuadrática y, a continuación, seleccione las raíces entre 0 y . La ecuación polinómica cuadrática en a resolver es con Por construcción y ; también aumenta en .Bx

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

A partir de esto, es fácil ver que si , la porción de parábola en la que estamos interesados ​​es una parte del lado derecho de la parábola, y la raíz a mantener es la más alta de las dos raíces, que es Por el contrario, si , la parábola está al revés y nos interesa su izquierda parte. La raíz a mantener es la más baja. Teniendo en cuenta el signo de , parece que esta es la misma raíz (es decir, la que tiene ) que en el primer caso.tanϕ>0tanϕ<0tanϕ

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

Aquí hay un código 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)))

histograma 1

Y con :ϕ<0

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

ingrese la descripción de la imagen aquí

Elvis
fuente
1
Cometí un error, porque establezco mi ángulo fuera de límites, lo entiendo. Pero su explicación de por qué debería evitar el uso de un solucionador numérico de todavía es confusa para mí. ¿Puedes intentar explicarlo más, por favor? Me encantaría conseguirlo. F(x)
Robert
@Robert Creo que su código funciona bien si el valor de es correcto. Sin embargo, le impide detectar posibles problemas (¿qué pasa si no hay una solución entre 0 y ? ¿O si ambas soluciones lo están? ¿O si no hay una solución real?). El trabajo adicional para evitar el uso del solucionador ya preparado vale la pena. BϕB
Elvis