Как реализовать точную перестановку в R? Опции

[ИНО]

Как посредством R случайным образом разделить ряд чисел на две группы заданных объемов с перебором всех возможных вариантов? Просто перемешивать исходный ряд и разделять в произвольной точке оказалось плохой идеей, так как порождаются в том числе и лишние варианты, когда состав обеих групп одинаков, но последовательность внутри них разная, например в итерации i в группе А оказываются 1, 2, 3, в группе Б - 4, 5, 6, а в итерации i+1 - соответственно 2, 1, 3 и 4, 5, 6. Если мы считаем статистику, инвариантную к порядку чисел, то эти паразитические псевдовыборки никак не сказываются на выводе, но сильно жрут вычислительные ресурсы. В общем, надо чтобы в каждой итерации группа А менялась с группой Б парой чисел.

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

[Игорь]

Напоминает критерий рандомизации Фишера-Питмана для связанных выборок.

[ИНО]

Не то, вариант для связанных не подходит, поскольку там первый член первой группы может поменяться только с первым членом второй группы и т. д. Мне же требуется, чтобы он мог поменяться также со вторым, третьим и т. д. Т. е. вариантов перестановки будет, больше чем для связанных выборок, но меньше чем при полном перемешивании. Причем это меньше выражается возможности или не возможности практической реализации начиная с определённого n. Так то полная рандомизация тоже сгодилось бы, но в итоге получится много дублей одного значения статистики. Например при разделении ряда из 6 разных чисел по двум равным группам выходит 72 дубля каждой величины межгрупповой разности средних. Если из каждого повторяющегося значения оставить только одно, то точное p-значения теста будет ровно таким же, как при сохранении всех, а вычислительные затраты сократятся в 72 раза! Сейчас большинство пакетов вычисляет полный перебор только для вектора длиной не более 7. Соответственно, при сбалансированном дизайне это будет по 3 в каждой группе. В таком случае независимо от величины различий минимальный достигаемый уровень значимости равен 0,1, что не может не печалить. Но точные тесты пакетах, например coin, обрабатывают и выборки побольше, притом шустро, думается их авторы секрет ограничения перестановки знают. Беда лишь в том, что эти тесты выдают только p, а мне требуются вся матрица перестановок.

В принципе, есть некоторые мысли как это запрограммировать с нуля: провести парные обмены как в случае связанных выборок, затем переставить значения только в второй группе и повторить, снова переставить только во второй группе и повторить, и т. д. Но есть подозрение что это будет громоздко и долго (вложенные циклы). Прежде, чем писать такое на R, хочется убедиться в отсутствии готовых более элегантных и шустрых решений.

[Игорь]

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

На больших выборках хорошо работает асимптотика.

В принципе, есть некоторые мысли как это запрограммировать с нуля: провести парные обмены как в случае связанных выборок, затем переставить значения только в второй группе и повторить, снова переставить только во второй группе и повторить, и т. д. Но есть подозрение что это будет громоздко и долго (вложенные циклы). Прежде, чем писать такое на R, хочется убедиться в отсутствии готовых более элегантных и шустрых решений.

Если бы Вы точнее поставили задачу и по возможности дали ссылку на источник, возникли бы более плодотворные идеи, как это реализовать. Например, перестановки - вроде бы понятная задача. А как получить сами перестановки, особенно для большой задачи, как Вы верно отметили, с множеством вложенных циклов? Пришлось сделать новый быстрый способ (хотя и требующий много памяти - но это сейчас не проблема) на основе матрицы планирования эксперимента (не опубликовано, исходники программы открыты).

Сообщение отредактировал Игорь - 30.05.2025 - 09:24

[ИНО]

Ссылку на источник чего? На идею, что при простой рандомизации объединенной последовательности появляются бесполезные варианта? Этот источник находится в моей голове и подтвержден результатами подсчетов.

По-моему я ясно все сформулировал, но попробую с другого ракурса. Задачу модно переформулировать как получение всех возможных перестановок чисел между двумя группами, но без бесполезных перестановок внутри каждой из групп. Например, если исходная группа А содержит числа 1, 2, а исходная группа Б - числа 10, 20, то вот полезные варианты перестановок:

1) А: 1, 10 и Б: 2, 20
2) А: 2, 10 и Б: 1, 20
3) А: 1, 20 и Б: 2, 10
4) А: 2, 20 и Б: 1, 10
5) А: 10, 20 и Б: 1, 2

А вот пример бесполезных (в случае, если уже имеются все вышеперечисленные) вариантов:

6) А: 2, 1 и Б: 20, 10 - для подсчета статистики, например, разности средних, это ровно то же , что и исходные группы.
7) А: 20, 2 и Б: 10, 1 - а это полностью эквивалентно варианту 4...

Все не буду перечислять, много их.

Так вот, обычный метод свободной рандомизации, описанный во множестве руководств, когда группы А и Б объединяют, объединенную последовательность тасуют всеми возможными способами и снова рвут на две группы посередине, порождает как полезные так и бесполезные перестановки. Причем бесполезных в разы больше. Уверен, что способ, дающий только полезные давным-давно придуман и реализован, но не умею искать.

[Игорь]

Ссылку на источник чего? На идею, что при простой рандомизации объединенной последовательности появляются бесполезные варианта?

При программировании критериев рандомизации я тоже "наткнулся" на эту проблему (по этой причине и "засветился" в данной теме) - слишком большой объем вычислений. Пришлось даже принудительно ограничивать возможности программы (не ждать же пользователю диалоговй программы, скажем, полчаса, пока расчет завершится). И подумал, нельзя ли исключить вычисление якобы бесполезных перестановок. Уже упоминал выше про "изобретенный" мной алгоритм (может быть, и изобретенный, т.к. ссылок не встречал, но не исключаю, что плохо искал, поэтому в кавычках) на основе матрицы планирования эксперимента 1-го порядка (мы ее использовали еще в 80-х годах для систематизации вычислений в одной специфической научно-технической области; есть много хороших источников, могу показать). Была идея использовать дробную реплику данной матрицы, что позволило бы существенно снизить объем вычислений. Оказалось, что в моих задачах все не совсем так. Это оказались не бесполезные перестановки - они есть и должны учитываться, даже если их вклад одинаков. Хотя в Вашей задаче, возможно, всё иначе. Поэтому и предложил показать не как бы только численную часть метода, а как это принято, начиная с постановки задачи, начиная со статистической модели. Но впрочем, это Ваша прерогатива.
Upd. То, что Вас интересует, похоже не на перестановки, а на сочетания.

Сообщение отредактировал Игорь - 31.05.2025 - 12:01

[ИНО]

Да, это именно сочетания и есть. Просто моя терминология в области комбинаторики оставляет желать много лучшего. Спасибо за подсказку термина, попробую поискать готовые решения с его помощью.

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

[ИНО]

Сделал! И даже пакетов дополнительных не понадобилось.

Код

recomb<-function(x, y)
{
  k<-max(c(length(x), length(y)))
  xy<-c(x, y)
  n<-length(xy)
  id<-1:n
  combmat<-combn(id, k)
  recb<-function(X, xy)
    {
    x<-xy[X]
    y<-xy[-X]
    if(length(x)==length(y)) return (data.frame(x=x, y=y)) else return(list(x, y))
    }
  result<-apply(combmat, 2, recb, xy=xy)
  return(result)
}

#Пример:

x<-c(1, 2)
y<-c(10, 20)
recomb(x, y)

Выдача (в случае равных длин обоих последовательностей):

Код

[[1]]
  x  y
1 1 10
2 2 20

[[2]]
   x  y
1  1  2
2 10 20

[[3]]
   x  y
1  1  2
2 20 10

[[4]]
   x  y
1  2  1
2 10 20

[[5]]
   x  y
1  2  1
2 20 10

[[6]]
   x y
1 10 1
2 20 2

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

Осталось разобраться, принесет ли мне сделанное пользу. Но это уже другая история.

[ИНО]

Итак, другая история: это работает в тестировании гипотез! Вот код из предыдущего поста, слегка допиленный для задачи проверки значимости разности средних:

Код

recombmeandifftest<-function(x, y)
{
  k<-max(c(length(x), length(y)))
  xy<-c(x, y)
  n<-length(xy)
  id<-1:n
  combmat<-combn(id, k)
  recmeandiff<-function(X, xy)
  {
    x<-xy[X]
    y<-xy[-X]
    meandiff<-abs(mean(x)-mean(y))
    return(meandiff)
  }
  diffobs<-abs(mean(x)-mean(y))
  diffsim<-apply(combmat[,2:ncol(combmat)], 2, recmeandiff, xy=xy)
  res<-as.numeric(round(diffsim, digits=6)>=round(diffobs, digits=6))
  p<-(1+sum(res))/(ncol(combmat))
  return(list(diffobs, p))
}

Проверим, на синтетических данных, немного большего размера, чем в прошлый раз, дабы не упереться в техническое препятствие (С)$100:

Код

x<-rnorm(10)
y<-rnorm(10, 3)
recombmeandifftest(x, y)

В выдаче первый элемент списка - наблюдаемый модуль разности средних (т. е. без перестановок), в торой - p-значение.

Сравним с тем что выдает критерий рандомизации компонент Фишера-Питмана, реaлизованный в пакете coin. Для этого придется преобразовать входные данные в табличную форму (другой их величество не принимают):

Код

xy<-data.frame(value=c(x, y), group=factor(c(rep("x", 10), rep("y", 10))))
oneway_test(value~group, data=xy, distribution="exact")

Я не смог наяндексить, что за статистику использует этот тест, и буду рад пояснению, что такое Z. С разностью средних не совпадает. Однако обратите внимание на практически идеальное совпадение p-значений! Напрашивается вывод: перестановки при полной рандомизации, которые я изначально посчитал лишними, таковыми и оказались. То есть, если объемы сравниваемых выборок равны между собой и равны n, то требуемое число перестановок, включая тождественную, - не вселяющий ужас (n*2)!, а всего лишь (n*2)!/(n!)^2, что гораздо веселее. С неравными выборками пока еще не разобрался.

[ИНО]

При неравных объемах выборок обнаружился баг - чувствительность к тому, что назначать х, а что у. Для правильного результата (соответствующего результатом oneway_test() из coin, требуется чтобы икс был длиннее игрека. Так и не понял в чем загвоздка но поваксил, тупо поменяв х и у местами в случае если второй длиннее. Принимайте обновленную версию функции:

Код

recombmeandifftest2<-function(x, y)
{
  if(length(y)>length(x))
     {
      x_<-y
      y_<-x}
     else
     {
      x_<-x
      y_<-y
     }
  k<-length(x_)
  xy<-c(x_, y_)
  n<-length(xy)
  id<-1:n
  combmat<-combn(id, k)
  recmeandiff<-function(X, xy)
  {
    x<-xy[X]
    y<-xy[-X]
    meandiff<-abs(mean(x)-mean(y))
    return(meandiff)
  }
  diffobs<-abs(mean(x_)-mean(y_))
  diffsim<-apply(combmat[,2:ncol(combmat)], 2, recmeandiff, xy=xy)
  res<-as.numeric(round(diffsim, digits=6)>=round(diffobs, digits=6))
  p<-(1+sum(res))/(ncol(combmat))
  return(list(diffobs, p))
}

Недостаток лишь один - думает медленно и жрет много памяти. Понятно, что пакет coin гораздо вычислительно эффективнее, чем моя самоделка, тем не менее не верю, что он генерирует все (n1+n2)! перестановок, скорее ограничивается все теми же уникальными сочетаниями, что и мой код. Потому как в случае с n1=15 и n2=10 на все про все у oneway_test() уходит лишь пара секунд. Для сравнения при вышеуказанных объемах выборок это 3268760 при сочетаниях и 15511210043330986055080688 при полной рандомизации. Очевидно, что второе вычислительно недостижимо, по крайней мере за те пару секунд Так что бесполезные варианты таки действительно бесполезны и, более того, ощутимо вредны!

Игорь, а в Вашем ПО количество перестановок получается какое?

Сообщение отредактировал ИНО - 1.06.2025 - 17:38