Вопрос по random, r – Генерация 3 случайных чисел, сумма которых равна 1 в R
Я надеюсь создать 3 (неотрицательных) квазислучайных числа, которые суммируются в одно и повторяются снова и снова.
В основном я пытаюсь разделить что-то на три случайные части на протяжении многих испытаний.
Пока я в курсе
a = runif(3,0,1)
Я думал, что я мог бы использовать1-a
как максимум в следующемrunif
, но это кажется грязным.
Но они, конечно, не сводятся к одному. Есть мысли, о, мудрые стекировщики?
это зависит от того, какое распределение вы хотите по номерам, но вот один из способов:
diff(c(0, sort(runif(2)), 1))
использованиеreplicate
чтобы получить столько наборов, сколько вы хотите:
> x <- replicate(5, diff(c(0, sort(runif(2)), 1)))
> x
[,1] [,2] [,3] [,4] [,5]
[1,] 0.66855903 0.01338052 0.3722026 0.4299087 0.67537181
[2,] 0.32130979 0.69666871 0.2670380 0.3359640 0.25860581
[3,] 0.01013117 0.28995078 0.3607594 0.2341273 0.06602238
> colSums(x)
[1] 1 1 1 1 1
a
а такжеb
тогда вы получили:
rand1 = min(a, b)
rand2 = abs(a - b)
rand3 = 1 - max(a, b)
которые прибавляют к 1 (или некоторому другому значению), тогда вам следует посмотреть наДирихле Дистрибьюшн.
Естьrdirichlet
функция вgtools
пакет и работаетRSiteSearch('Dirichlet')
поднимает немало хитов, которые могут легко привести вас к инструментам для этого (и нетрудно написать код вручную для простых дистрибутивов Дирихле).
чем может показаться на первый взгляд. Посмотрев на следующее, вы можете тщательно продумать процесс, который вы используете для представления этих чисел:
## My initial idea (and commenter Anders Gustafsson's):
## Sample 3 random numbers from [0,1], sum them, and normalize
jobFun <- function(n) {
m <- matrix(runif(3*n,0,1), ncol=3)
m<- sweep(m, 1, rowSums(m), FUN="/")
m
}
## Andrie's solution. Sample 1 number from [0,1], then break upper
## interval in two. (aka "Broken stick" distribution).
andFun <- function(n){
x1 <- runif(n)
x2 <- runif(n)*(1-x1)
matrix(c(x1, x2, 1-(x1+x2)), ncol=3)
}
## ddzialak's solution (vectorized by me)
ddzFun <- function(n) {
a <- runif(n, 0, 1)
b <- runif(n, 0, 1)
rand1 = pmin(a, b)
rand2 = abs(a - b)
rand3 = 1 - pmax(a, b)
cbind(rand1, rand2, rand3)
}
## Simulate 10k triplets using each of the functions above
JOB <- jobFun(10000)
AND <- andFun(10000)
DDZ <- ddzFun(10000)
## Plot the distributions of values
par(mfcol=c(2,2))
hist(JOB, main="JOB")
hist(AND, main="AND")
hist(DDZ, main="DDZ")
ьшой тест из трех предложенных основных алгоритмов и того, какие средние значения они дадут для сгенерированных чисел.
choose_one_and_divide_rest
means: [ 0.49999212 0.24982403 0.25018384]
standard deviations: [ 0.28849948 0.22032758 0.22049302]
time needed to fill array of size 1000000 was 26.874945879 seconds
choose_two_points_and_use_intervals
means: [ 0.33301421 0.33392816 0.33305763]
standard deviations: [ 0.23565652 0.23579615 0.23554689]
time needed to fill array of size 1000000 was 28.8600130081 seconds
choose_three_and_normalize
means: [ 0.33334531 0.33336692 0.33328777]
standard deviations: [ 0.17964206 0.17974085 0.17968462]
time needed to fill array of size 1000000 was 27.4301018715 seconds
Измерения времени должны выполняться с большой долей соли, так как они могут в большей степени зависеть от управления памятью Python, чем от самого алгоритма. Мне лень делать это правильно сtimeit
, Я сделал это на 1GHz Atom, так что это объясняет, почему это заняло так много времени.
Во всяком случае, choose_one_and_divide_rest - это алгоритм, предложенный Андри, и сам плакат с вопросом (А ТАКЖЕ): вы выбираете одно значение a в [0,1], затем одно в [a, 1] и затем смотрите, что у вас осталось. Это добавляет к одному, но это об этом, первое деление в два раза больше, чем два других. Можно было бы догадаться, как много ...
choose_two_points_and_use_intervals - это принятый ответ от ddzialak (ДДЗ). Он принимает две точки в интервале [0,1] и использует размер трех подинтервалов, созданных этими точками, в качестве трех чисел. Работает как шарм и средства все 1/3.
Choose_three_and_normalize - это решение, разработанное Андерсом Густафссоном и Джошем О'Брайеном (РАБОТА). Это просто g, генерирует три числа в [0,1] и нормализует их обратно к сумме 1. Работает так же хорошо и на удивление немного быстрее в моей реализации Python. Дисперсия немного ниже, чем для второго решения.
Там у вас есть это. Не знаю, какому бета-распределению соответствуют эти решения или какой набор параметров в соответствующей статье я упоминал в комментарии, но, возможно, кто-то другой может это выяснить.