Ideal () в пакете R pscl не дает повторяющихся результатов

Как подозревают OP и @SchaunW, проблема заключается в коде C. «Немного» копания выявило довольно тонкую проблему (см. источник код, хотя и не самая новая версия):

Вся выборка в Ideal.c появляется в части начало итераций, т. е. там, где используются функции updatex, updatey и другие. Однако проблема с одним из аргументов этих функций — матрицей ok (иронично, правда?). Он используется updatex и updateb и только позициями, где важны ok == 1 (в crosscheck, crosscheckx).

До этого некоторым значениям ok присваивается значение 1 в check(y,ok,n,m).

Однако в самом начале начальные значения ok обозначаются

ok = imatrix(n,m);

который выделяет целочисленную матрицу (см. util.c для imatrix). Проблема в том, что тогда ok содержит разные числа, т.е. не только нули, а иногда и единицы. Кажется, что они не связаны с состоянием RNG R, что объясняет поведение, отмеченное @SchaunW: all.equal(run1@result$xbar, run2@result$xbar) возвращает TRUE, если !any(ok == 1), и наоборот. Кроме того, разное количество единиц объясняет разное lastState.

Я не эксперт в C, и я не уверен, есть ли логическая ошибка в коде или функция imatrix должна быть исправлена, но простым решением может быть заполнение ok нулями сразу после инициализации:

ok = imatrix(n,m);
for(a=0; a<n; a++) {
    for(aa=0; aa<m; aa++) {
      ok[a][aa] = 0;
    }
}

Наконец, есть также исправление, не включающее изменение кода C (хотя оно может не подходить для ваших приложений). Функции crossxyi, crossxyj используются вместо crosscheck, crosscheckx (плохих), когда impute = TRUE вместо ideal.

ИЗМЕНИТЬ

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

Проблема, кажется, заключается в коде C. Когда я открываю функцию ideal и запускаю ее построчно, all.equal возвращает TRUE для каждого ввода в этой строке кода:

output <- .C("IDEAL", PACKAGE = .package.Name, as.integer(n), 
      as.integer(m), as.integer(d), as.double(yToC), as.integer(maxiter), 
      as.integer(thin), as.integer(impute), as.integer(mda), 
      as.double(xp), as.double(xpv), as.double(bp), as.double(bpv), 
      as.double(xstart), as.double(bstart), xoutput = as.double(rep(0, 
        n * d * numrec)), boutput = as.double(0), as.integer(burnin), 
      as.integer(usefile), as.integer(store.item), as.character(file), 
      as.integer(verbose))

Однако, когда я запускаю приведенный выше код несколько раз, output$xoutput каждый раз возвращает несколько разные результаты, даже если я вызываю set.seed(42) непосредственно перед каждым запуском.

sessionInfo()
R version 2.15.2 (2012-10-26)
Platform: x86_64-apple-darwin9.8.0/x86_64 (64-bit)

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] grid      splines   stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] SoDA_1.0-5       pscl_1.04.4      vcd_1.2-13       colorspace_1.2-0 gam_1.06.2       coda_0.16-1      lattice_0.20-10  mvtnorm_0.9-9994
[9] MASS_7.3-22     

loaded via a namespace (and not attached):
[1] tools_2.15.2

ИСХОДНЫЙ ОТВЕТ

Функция ideal имеет аргумент startvals. Значение по умолчанию для этого аргумента — «eigen». Чтобы ваш вызов set.seed имел эффект, вам нужно изменить этот аргумент на «случайный». Вот что вы уже пробовали:

run1 <- simulationResult(
   ideal(s109, 
     normalize=TRUE,
     maxiter = 500,
     thin = 10,
     burnin = 0,
     startvals = "eigen"),
   seed = 42)

run2 <- simulationResult(
   ideal(s109, 
     normalize=TRUE,
     maxiter = 500,
     thin = 10,
     burnin = 0,
     startvals = "eigen"),
   seed = 42)

all.equal(run1@firstState, run2@firstState)
[1] TRUE

all.equal(run1@result$xbar, run2@result$xbar)
[1] "Mean relative difference: 0.01832379"

А вот то же самое с startvals, установленным на «случайный»:

run1 <- simulationResult(
   ideal(s109, 
     normalize=TRUE,
     maxiter = 500,
     thin = 10,
     burnin = 0,
     startvals = "random"),
   seed = 42)

run2 <- simulationResult(
   ideal(s109, 
     normalize=TRUE,
     maxiter = 500,
     thin = 10,
     burnin = 0,
     startvals = "random"),
   seed = 42)

all.equal(run1@firstState, run2@firstState)
[1] TRUE    

all.equal(run1@result$xbar, run2@result$xbar)
[1] TRUE

Насколько я вижу, необходимость установки startvals в "random" для получения воспроизводимых результатов четко не указана в документации пакета. Мне пришлось поиграть с ним некоторое время, прежде чем я понял это.

Это модель MCMC, поэтому она обязательно использует генерацию случайных чисел. Чтобы получить воспроизводимые результаты, вам нужно начать анализ, установив «начальное число» для генератора случайных чисел. Таким образом, каждый раз, когда вы строите модель, она использует одни и те же «случайные» числа (пока вы сбрасываете начальное число каждый раз, когда строите модель. Используйте функцию set.seed() и просто вводите в нее произвольное значение, например 1234.

Я не знаком с этим пакетом, и похоже, что вы, возможно, уже устанавливаете начальное число для генерации случайных чисел в своем вызове функции с помощью seed=42, но я бы рекомендовал установить его явно с помощью set.seed() в любом случае. Затем ваш код становится:

set.seed(1234)
run1 <- simulationResult(
  ideal(s109, 
    normalize=TRUE,
    maxiter = 500,
    thin = 10,
    burnin = 0),
  seed = 42)

set.seed(1234)
run2 <- simulationResult(
  ideal(s109, 
    normalize=TRUE,
    maxiter = 500,
    thin = 10,
    burnin = 0),
  seed = 42)