Генерация случайных чисел из «наклонного равномерного распределения» из математической теории

9

Для каких-то целей мне нужно генерировать случайные числа (данные) из распределения "наклонной формы". «Наклон» этого распределения может изменяться в некотором разумном интервале, и тогда мое распределение должно измениться с равномерного на треугольное в зависимости от наклона. Вот мой вывод:

введите описание изображения здесь

Давайте сделаем это просто и сгенерируем данные от до (синий, красный - равномерное распределение). Чтобы получить функцию плотности вероятности синей линии, мне нужно просто уравнение этой линии. Таким образом:0B

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

и с тех пор (фото):

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

У нас есть это:

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

Поскольку - это PDF, CDF равен:f(x)

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

Теперь давайте сделаем генератор данных. Идея состоит в том, что если я исправлю , случайные числа могут быть вычислены, если я получу числа из из равномерного распределения, как описано здесь . Таким образом, если мне нужно 100 случайных чисел из моего распределения с фиксированным , то для любого из равномерного распределения существует из «наклонного распределения», и можно вычислить как:φ,Bx(0,1)φ,Bti(0,1)xix

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

Из этой теории я сделал код на Python, который выглядит следующим образом:

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

Но сгенерированные числа rand_numbочень близки к нулю или к B (который я установил как 25). Разницы нет, когда я генерирую 100 чисел, все они близки к 25 или все близки к нулю. В один заход:

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]

Так что в моем коде должно быть что-то не так. Может ли кто-нибудь помочь мне с моим выводом или кодом? Я без ума от этого сейчас, я не вижу никакой ошибки. Я полагаю, код R даст мне аналогичные результаты.

Роберт
источник
2
Если вам нужно только генерировать случайные числа, вам вообще не нужно определять распределение. Просто бросьте дротики на свое изображение и сохраните их x-координаты, но когда дротик приземлится в левом треугольнике, обозначенном « », измените его x-координату с на . Например, присвойте любые значения и (реальный параметр, который при заданных значениях от до создает ваши распределения) и задайте количество необходимых вам случайных значений. Вот код:ϕxBxBtheta11nRx<-runif(n,-1,1);x<-(ifelse(runif(n,-1,1)>theta*x,-x,x)+1)*(B/2)
whuber

Ответы:

9

Ваш вывод в порядке. Обратите внимание, что для получения положительной плотности , вы должны ограничить В вашем коде поэтому вы должны взять между , вот где ваш код не работает.(0,B)

B2tanϕ<2.
B=25ϕ±tan12625

Вы можете (и должны) избегать использования квадратичного решателя, а затем выбрать корни между 0 и . Квадратичное полиномиальное уравнение в должно быть решено, имеет вид с По построению и ; также увеличивается на .Bx

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

Из этого легко видеть, что если , то часть параболы, в которой мы заинтересованы, является частью правой стороны параболы, а корень, который нужно сохранить, является самым высоким из двух корней, что is Наоборот, если , парабола перевернута, и мы заинтересованы в ее левой часть. Корень, который нужно сохранить - самый низкий. Принимая во внимание знак кажется, что это тот же корень (т. у которого ), чем в первом случае.tanϕ>0

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

Вот немного кода 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)))

гистограмма 1

И с :ϕ<0

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

введите описание изображения здесь

Элвис
источник
1
Я сделал ошибку, потому что я установил свой угол за пределы, я понял. Но ваше объяснение, почему я должен избегать использования числового решателя для меня все еще туманно. Можете ли вы попытаться объяснить это больше, пожалуйста? Я люблю, чтобы получить это. F(x)
Роберт
@ Роберт Я думаю, твой код работает хорошо, если значение верное. Тем не менее, он не дает вам уловить потенциальные проблемы (что, если между 0 и нет решения ? Или если оба решения есть? Или нет реального решения?). Дополнительная работа, чтобы избежать использования готового решателя, стоит того. ϕB
Элвис