Назначение определенного количества значений, основанных на распределении вероятностей (в R) - PullRequest
7 голосов
/ 04 августа 2011

Привет и заранее спасибо за помощь!

Я пытаюсь сгенерировать вектор с определенным количеством значений, которые назначаются в соответствии с распределением вероятностей. Например, я хочу вектор длиной 31, содержащий 26 нулей и 5 единиц. (Общая сумма вектора должна всегда быть пятью.) Однако, расположение их важно. И чтобы определить, какие значения должны быть равны одному, а какие - нулю, у меня есть вектор вероятностей (длина 31), который выглядит так:

probs<-c(0.01,0.02,0.01,0.02,0.01,0.01,0.01,0.04,0.01,0.01,0.12,0.01,0.02,0.01,
0.14,0.06,0.01,0.01,0.01,0.01,0.01,0.14,0.01,0.07,0.01,0.01,0.04,0.08,0.01,0.02,0.01)

Я могу выбрать значения в соответствии с этим распределением и получить вектор длины 31, используя rbinom, но я не могу выбрать ровно пять значений.

Inv=rbinom(length(probs),1,probs)
Inv
[1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0

Есть идеи?

Еще раз спасибо!

Ответы [ 3 ]

10 голосов
/ 04 августа 2011

Как насчет использования взвешенного sample.int для выбора местоположений?

Inv<-integer(31)
Inv[sample.int(31,5,prob=probs)]<-1
Inv
[1] 0 0 0 1 0 1 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0
7 голосов
/ 04 августа 2011

Чейз дает отличный ответ и упоминает проблему повторения while().Одна из проблем, связанных с побегом while(), заключается в том, что если вы делаете это одно испытание за раз, и потребуется много, скажем t , испытаний, чтобы найти то, которое соответствует целевому числу 1 s, вы получаете накладные расходы на t вызовов основной функции, в данном случае rbinom().

Однако есть выход, потому что rbinom(), как и всеэти (псевдо) генераторы случайных чисел в R векторизованы, мы можем генерировать m испытаний одновременно и проверять эти m испытаний на соответствие требованиям 5 1 с.Если ничего не найдено, мы неоднократно проводим m испытаний, пока не найдем то, которое соответствует.Эта идея реализована в функции foo() ниже.Аргумент chunkSize равен m , количество попыток, которые нужно провести за раз.Я также воспользовался возможностью, чтобы позволить функции найти больше, чем одно конформное испытание;Аргумент n контролирует количество возвращаемых конформных испытаний.

foo <- function(probs, target, n = 1, chunkSize = 100) {
    len <- length(probs)
    out <- matrix(ncol = len, nrow = 0) ## return object
    ## draw chunkSize trials
    trial <- matrix(rbinom(len * chunkSize, 1, probs),
                    ncol = len, byrow = TRUE)
    rs <- rowSums(trial)  ## How manys `1`s
    ok <- which(rs == 5L) ## which meet the `target`
    found <- length(ok)   ## how many meet the target
    if(found > 0)         ## if we found some, add them to out
        out <- rbind(out,
                     trial[ok, , drop = FALSE][seq_len(min(n,found)), , 
                                               drop = FALSE])
    ## if we haven't found enough, repeat the whole thing until we do
    while(found < n) {
        trial <- matrix(rbinom(len * chunkSize, 1, probs),
                            ncol = len, byrow = TRUE)
        rs <- rowSums(trial)
        ok <- which(rs == 5L)
        New <- length(ok)
        if(New > 0) {
            found <- found + New
            out <- rbind(out, trial[ok, , drop = FALSE][seq_len(min(n, New)), , 
                                                        drop = FALSE])
        }
    }
    if(n == 1L)           ## comment this, and
        out <- drop(out)  ## this if you don't want dimension dropping
    out
}

Он работает так:

> set.seed(1)
> foo(probs, target = 5)
 [1] 1 0 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0
[31] 0
> foo(probs, target = 5, n = 2)
     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
[1,]    0    0    0    0    0    0    0    0    0     0     0
[2,]    0    0    0    0    0    0    0    0    0     0     1
     [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21]
[1,]     0     0     0     1     1     0     0     0     0     0
[2,]     0     1     0     0     1     0     0     0     0     0
     [,22] [,23] [,24] [,25] [,26] [,27] [,28] [,29] [,30] [,31]
[1,]     1     0     1     0     0     0     1     0     0     0
[2,]     1     0     1     0     0     0     0     0     0     0

Обратите внимание, что я отбрасываю пустое измерение в случае, когда n == 1.Прокомментируйте последний фрагмент кода if, если вы не хотите использовать эту функцию.

Вам необходимо сбалансировать размер chunkSize с вычислительной нагрузкой, связанной с проверкой такого количества испытаний одновременно.Если требование (здесь 5 1 с) очень маловероятно, увеличьте chunkSize, чтобы уменьшить количество вызовов до rbinom().Если требование является вероятным, есть небольшие испытания для рисования точек и большие chunkSize за один раз, если вы хотите только одно или два, так как вы должны оценивать каждую пробную игру.

5 голосов
/ 04 августа 2011

Я думаю, что вы хотите выполнить повторную выборку из биномиального распределения с заданным набором вероятностей, пока не достигнете целевого значения 5, это верно? Если так, то я думаю, что это делает то, что вы хотите. Цикл while может использоваться для итерации до тех пор, пока условие не будет выполнено. Если вы вводите очень нереалистичные вероятности и целевые значения, я думаю, это может превратиться в убегающую функцию, так что считайте себя предупрежденным:)

FOO <- function(probs, target) {
  out <- rbinom(length(probs), 1, probs)

  while (sum(out) != target) {

    out <- rbinom(length(probs), 1, probs)
  }
  return(out)
}

FOO (пробники, цель = 5)

> FOO(probs, target = 5)  
 [1] 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 1 0 0 1 0 0 0 0 1 0
...