Алгоритм отбора проб без замены?

Question

Я пытаюсь проверить вероятность того, что определенная кластеризация данных произошла случайно. Надежным способом сделать это является моделирование по методу Монте-Карло, в котором ассоциации между данными и группами случайным образом переназначаются большое количество раз (например, 10 000), а метрика кластеризации используется для сравнения фактических данных с моделированием для определения значение.

У меня большая часть работы с указателями, отображающими группировку на элементы данных, поэтому я планирую случайным образом переназначить указатели на данные. ВОПРОС: что такое быстрый способ выборки без замены, чтобы каждый указатель случайным образом переназначался в повторяющихся наборах данных?

Например (эти данные являются лишь упрощенным примером):

Данные (n = 12 значений) - Группа A: 0,1, 0,2, 0,4 / Группа B: 0,5, 0,6, 0,8 / Группа C: 0,4, 0,5 / Группа D: 0,2, 0,2, 0,3, 0,5

Для каждого дублированного набора данных у меня будут одинаковые размеры кластеров (A = 3, B = 3, C = 2, D = 4) и значения данных, но я переназначу значения кластерам.

Для этого я мог бы генерировать случайные числа в диапазоне 1-12, назначать первый элемент группы A, затем генерировать случайные числа в диапазоне 1-11 и назначать второй элемент в группе A и так далее. Переназначение указателя происходит быстро, и я предварительно выделю все структуры данных, но выборка без замены кажется проблемой, которую можно было решить много раз раньше.

Предпочтителен логический или псевдокод.

Answer 1

Вот некоторый код для выборки без замены, основанный на Алгоритме 3.4.2S книги Кнута «Полу-цифровые алгоритмы».

void SampleWithoutReplacement
(
    int populationSize,    // size of set sampling from
    int sampleSize,        // size of each sample
    vector<int> & samples  // output, zero-offset indicies to selected items
)
{
    // Use Knuth's variable names
    int& n = sampleSize;
    int& N = populationSize;

    int t = 0; // total input records dealt with
    int m = 0; // number of items selected so far
    double u;

    while (m < n)
    {
        u = GetUniform(); // call a uniform(0,1) random number generator

        if ( (N - t)*u >= n - m )
        {
            t++;
        }
        else
        {
            samples[m] = t;
            t++; m++;
        }
    }
}

Существует более эффективный, но более сложный метод Джеффри Скотта Виттера в "Эффективном алгоритме последовательной случайной выборки", транзакции ACM по математическому программному обеспечению, 13 (1), март 1987 г., 58-67.

Answer 2

Рабочий код C ++, основанный на ответе Джона Д. Кука .

#include <random>
#include <vector>

double GetUniform()
{
    static std::default_random_engine re;
    static std::uniform_real_distribution<double> Dist(0,1);
    return Dist(re);
}

// John D. Cook, https://stackoverflow.com/a/311716/15485
void SampleWithoutReplacement
(
    int populationSize,    // size of set sampling from
    int sampleSize,        // size of each sample
    std::vector<int> & samples  // output, zero-offset indicies to selected items
)
{
    // Use Knuth's variable names
    int& n = sampleSize;
    int& N = populationSize;

    int t = 0; // total input records dealt with
    int m = 0; // number of items selected so far
    double u;

    while (m < n)
    {
        u = GetUniform(); // call a uniform(0,1) random number generator

        if ( (N - t)*u >= n - m )
        {
            t++;
        }
        else
        {
            samples[m] = t;
            t++; m++;
        }
    }
}

#include <iostream>
int main(int,char**)
{
  const size_t sz = 10;
  std::vector< int > samples(sz);
  SampleWithoutReplacement(10*sz,sz,samples);
  for (size_t i = 0; i < sz; i++ ) {
    std::cout << samples[i] << "\t";
  }

  return 0;
}

Answer 3

См. Мой ответ на этот вопрос Уникальные (неповторяющиеся) случайные числа в O (1)? . Та же логика должна выполнить то, что вы хотите сделать.

Answer 4

Вдохновленный @ ответом Джона Д. Кука , я написал реализацию в Nim. Сначала мне было трудно понять, как это работает, поэтому я подробно прокомментировал, включая пример. Может быть, это помогает понять идею. Также я немного изменил имена переменных.

iterator uniqueRandomValuesBelow*(N, M: int) =
  ## Returns a total of M unique random values i with 0 <= i < N
  ## These indices can be used to construct e.g. a random sample without replacement
  assert(M <= N)

  var t = 0 # total input records dealt with
  var m = 0 # number of items selected so far

  while (m < M):
    let u = random(1.0) # call a uniform(0,1) random number generator

    # meaning of the following terms:
    # (N - t) is the total number of remaining draws left (initially just N)
    # (M - m) is the number how many of these remaining draw must be positive (initially just M)
    # => Probability for next draw = (M-m) / (N-t)
    #    i.e.: (required positive draws left) / (total draw left)
    #
    # This is implemented by the inequality expression below:
    # - the larger (M-m), the larger the probability of a positive draw
    # - for (N-t) == (M-m), the term on the left is always smaller => we will draw 100%
    # - for (N-t) >> (M-m), we must get a very small u
    #
    # example: (N-t) = 7, (M-m) = 5
    # => we draw the next with prob 5/7
    #    lets assume the draw fails
    # => t += 1 => (N-t) = 6
    # => we draw the next with prob 5/6
    #    lets assume the draw succeeds
    # => t += 1, m += 1 => (N-t) = 5, (M-m) = 4
    # => we draw the next with prob 4/5
    #    lets assume the draw fails
    # => t += 1 => (N-t) = 4
    # => we draw the next with prob 4/4, i.e.,
    #    we will draw with certainty from now on
    #    (in the next steps we get prob 3/3, 2/2, ...)
    if (N - t)*u >= (M - m).toFloat: # this is essentially a draw with P = (M-m) / (N-t)
      # no draw -- happens mainly for (N-t) >> (M-m) and/or high u
      t += 1
    else:
      # draw t -- happens when (M-m) gets large and/or low u
      yield t # this is where we output an index, can be used to sample
      t += 1
      m += 1

# example use
for i in uniqueRandomValuesBelow(100, 5):
  echo i

Answer 5

Когда размер популяции намного превышает размер выборки, вышеприведенные алгоритмы становятся неэффективными, поскольку имеют сложность O ( n ), n Численность населения.

Когда я был студентом, я написал несколько алгоритмов равномерной выборки без замены, которые имеют среднюю сложность O ( s log s ), где s - размер выборки. Вот код для алгоритма двоичного дерева, со средней сложностью O ( s log s ), в R:

# The Tree growing algorithm for uniform sampling without replacement
# by Pavel Ruzankin 
quicksample = function (n,size)
# n - the number of items to choose from
# size - the sample size
{
  s=as.integer(size)
  if (s>n) {
    stop("Sample size is greater than the number of items to choose from")
  }
  # upv=integer(s) #level up edge is pointing to
  leftv=integer(s) #left edge is poiting to; must be filled with zeros
  rightv=integer(s) #right edge is pointig to; must be filled with zeros
  samp=integer(s) #the sample
  ordn=integer(s) #relative ordinal number

  ordn[1L]=1L #initial value for the root vertex
  samp[1L]=sample(n,1L) 
  if (s > 1L) for (j in 2L:s) {
    curn=sample(n-j+1L,1L) #current number sampled
    curordn=0L #currend ordinal number
    v=1L #current vertice
    from=1L #how have come here: 0 - by left edge, 1 - by right edge
    repeat {
      curordn=curordn+ordn[v]
      if (curn+curordn>samp[v]) { #going down by the right edge
        if (from == 0L) {
          ordn[v]=ordn[v]-1L
        }
        if (rightv[v]!=0L) {
          v=rightv[v]
          from=1L
        } else { #creating a new vertex
          samp[j]=curn+curordn
          ordn[j]=1L
          # upv[j]=v
          rightv[v]=j
          break
        }
      } else { #going down by the left edge
        if (from==1L) {
          ordn[v]=ordn[v]+1L
        }
        if (leftv[v]!=0L) {
          v=leftv[v]
          from=0L
        } else { #creating a new vertex
          samp[j]=curn+curordn-1L
          ordn[j]=-1L
          # upv[j]=v
          leftv[v]=j
          break
        }
      }
    }
  }
  return(samp)  
}

Сложность этого алгоритма обсуждается в: Rouzankin, P. S .; Войтишек А. В. О стоимости алгоритмов случайного выбора. Применение методов Монте-Карло 5 (1999), № 1, 39-54. http://dx.doi.org/10.1515/mcma.1999.5.1.39

Если вы найдете алгоритм полезным, сделайте ссылку.

Смотри также: П. Гупта, Г. П. Бхаттачарджи. (1984) Эффективный алгоритм случайной выборки без замены. Международный журнал по компьютерной математике 16: 4, страницы 201-209. DOI: 10.1080 / 00207168408803438

Теухола, Дж. И Невалайнен, О. 1982. Два эффективных алгоритма для случайной выборки без замены. / IJCM /, 11 (2): 127–140. DOI: 10.1080 / 00207168208803304

В последней статье авторы используют хеш-таблицы и утверждают, что их алгоритмы имеют сложность O ( s ). Есть еще один быстрый алгоритм хеш-таблицы, который скоро будет реализован в pqR (довольно быстрый R): https://stat.ethz.ch/pipermail/r-devel/2017-October/075012.html

Answer 6

Другой алгоритм отбора проб без замены описан здесь .

Это похоже на то, что описал Джон Д. Кук в своем ответе, а также от Кнута, но у него другая гипотеза: размер популяции неизвестен, но выборка может уместиться в памяти. Этот называется «алгоритм Кнута S».

Цитирование статьи с розеттакодом:

1. Выберите первые n элементов в качестве образца по мере их появления;
2. Для i-го предмета, где i> n, есть случайный шанс n / i его сохранить. Если этот шанс не получится, образец останется прежним. Если нет, случайным образом (1 / n) заменить один из ранее выбранных n элементы образца.
3. Повторите # 2 для всех последующих пунктов.