Вопрос по random, r – Генерация 3 случайных чисел, сумма которых равна 1 в R

11

Я надеюсь создать 3 (неотрицательных) квазислучайных числа, которые суммируются в одно и повторяются снова и снова.

В основном я пытаюсь разделить что-то на три случайные части на протяжении многих испытаний.

Пока я в курсе

a = runif(3,0,1)

Я думал, что я мог бы использовать1-a как максимум в следующемrunif, но это кажется грязным.

Но они, конечно, не сводятся к одному. Есть мысли, о, мудрые стекировщики?

Как насчет генерации 2 случайных чисел a и b? Тогда a + b + c = 1 = & gt; с = 1 - (а + б) Frank Schmitt
Кажется, что ваша проблема называется «случайным делением интервала». Это классическая проблема в статистике, и поиск в Google дает удивительно маленькие красивые картинки ... Christian
а если a и b сумма больше 1? mmann1123
какое распределение вы хотите по интервалам? Например, должна ли быть равномерно выбрана длина первого интервала по (0,1)? ALiX
Есть ли возможность перенормировать случайные числа после генерации? Anders Gustafsson

Ваш Ответ

6   ответов
1

ения, а затем разделил бы их по сумме. Код как ниже.

n <- 3
x <- runif(3, 0, 1)
y <- x/sum(x)
sum(y)== 1

n может быть любым числом, которое вам нравится.

5

это зависит от того, какое распределение вы хотите по номерам, но вот один из способов:

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
10

a а такжеb тогда вы получили:

rand1 = min(a, b)
rand2 = abs(a - b)
rand3 = 1 - max(a, b)
@ пользователь так a = 0,85, b = 0,99, тогда вы получили числа: 0,85, 0,14, 0,01 (для меня это очень хорошие 3 случайных числа из 0..1)
Получающееся распределение, кажется, не совсем тривиально:jstor.org/discover/10.2307/… и более поздняя статья, которая может быть свободно доступнаdoc.utwente.nl/70657/1/Sleutel67random.pdf
Кроме того, вы должны повторить генерацию второго числа, если a == b ... (должно быть ОЧЕНЬ редкий случай)
8

которые прибавляют к 1 (или некоторому другому значению), тогда вам следует посмотреть наДирихле Дистрибьюшн.

Естьrdirichlet функция вgtools пакет и работаетRSiteSearch('Dirichlet') поднимает немало хитов, которые могут легко привести вас к инструментам для этого (и нетрудно написать код вручную для простых дистрибутивов Дирихле).

13

чем может показаться на первый взгляд. Посмотрев на следующее, вы можете тщательно продумать процесс, который вы используете для представления этих чисел:

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

enter image description here

Хорошо, я думал о графике результатов, но вы уже сделали это. Интересно видеть, что, очевидно, ни одно из решений на самом деле не делает то, что хотелось бы интуитивно. Интересно также, что на этих графиках вы действительно не видите, что DDZ делает правильные вещи в соответствии со средствами, в то время как AND даже не делает этого.
2

ьшой тест из трех предложенных основных алгоритмов и того, какие средние значения они дадут для сгенерированных чисел.

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. Дисперсия немного ниже, чем для второго решения.

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

Похожие вопросы