Для каких-то целей мне нужно генерировать случайные числа (данные) из распределения "наклонной формы". «Наклон» этого распределения может изменяться в некотором разумном интервале, и тогда мое распределение должно измениться с равномерного на треугольное в зависимости от наклона. Вот мой вывод:
Давайте сделаем это просто и сгенерируем данные от до (синий, красный - равномерное распределение). Чтобы получить функцию плотности вероятности синей линии, мне нужно просто уравнение этой линии. Таким образом:
и с тех пор (фото):
У нас есть это:
Поскольку - это PDF, CDF равен:
Теперь давайте сделаем генератор данных. Идея состоит в том, что если я исправлю , случайные числа могут быть вычислены, если я получу числа из из равномерного распределения, как описано здесь . Таким образом, если мне нужно 100 случайных чисел из моего распределения с фиксированным , то для любого из равномерного распределения существует из «наклонного распределения», и можно вычислить как:
Из этой теории я сделал код на 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 даст мне аналогичные результаты.
B
theta
n
R
x<-runif(n,-1,1);x<-(ifelse(runif(n,-1,1)>theta*x,-x,x)+1)*(B/2)
Ответы:
Ваш вывод в порядке. Обратите внимание, что для получения положительной плотности , вы должны ограничить В вашем коде поэтому вы должны взять между , вот где ваш код не работает.(0,B)
Вы можете (и должны) избегать использования квадратичного решателя, а затем выбрать корни между 0 и . Квадратичное полиномиальное уравнение в должно быть решено, имеет вид с По построению и ; также увеличивается на .B x
Из этого легко видеть, что если , то часть параболы, в которой мы заинтересованы, является частью правой стороны параболы, а корень, который нужно сохранить, является самым высоким из двух корней, что is Наоборот, если , парабола перевернута, и мы заинтересованы в ее левой часть. Корень, который нужно сохранить - самый низкий. Принимая во внимание знак кажется, что это тот же корень (т. у которого ), чем в первом случае.tanϕ>0
Вот немного кода R
И с :ϕ<0
источник