ChoiceModelR, none=TRUE, ошибка в betadraw[хорошо, ] = новаябета[хорошо, ]

Я использую ChoiceModelR для иерархического полиномиального логита. Я хочу получить оценки полезности внешнего товара (который соответствует нормальному распределению). У внешнего товара нет ковариатов, как у внутренних товаров, т.е. у него не может быть фиктивной цены или марки - , поэтому я устанавливаю list(none=TRUE) и не добавляю это отсутствие выбора к данным X (как описано в документации ChoiceModelR), а только к данным y (выбор).

Итерации запускаются нормально, потом в какой-то момент останавливается и говорит

"Error in betadraw[good, ] = newbeta[good, ] :   NAs are not allowed in subscripted assignments". 

Вероятно, это происходит из-за того, что в строке 388 функции "choicemodelr" индекс "хорошо" равен NA.

Я просмотрел несколько вопросов о selectionmodelr (это,это,this), а также о NA в нижнем индексе (это, это), но я предполагаю, что моя проблема специфична для этой функции в том смысле, что, вероятно, некоторые входные данные в итерации становятся настолько большими/маленькими, что «хорошо» превратится в NA.

Ниже приведен очень простой пример. Я генерирую данные с 3 продуктами с различными атрибутами. В половине периодов продукт 3 не предлагается. 2000 потребителей имеют предпочтения, распределенные нормально, по трем атрибутам (и предпочтение внешнего товара). Ошибка логита добавлена ​​для соответствия модели. Внешний товар индексируется как продукт 4 (и когда в наборе выбора было 3 и 2 продукта).

Как я мог избежать ошибки NA? Я что-то не так делаю, или это общий баг в функции?

Я также искал в Интернете примеры установки опции none=TRUE, но не нашел ни одного воспроизводимого. Возможно, этот вариант является только проблематичным, так как нет проблем с восстановлением истинных параметров, если я устанавливаю none=FALSE и не позволяю клиентам выбирать внешний вариант.

Таким образом, код, который приводит к ошибке NA, следующий:

library("ChoiceModelR")
library("MASS")

set.seed(36)

# Set demand pars
beta_mu = c(-3,4,1)
beta_sigma = diag(c(1,1,1))
alfa_mu = 5  #outside good mean utility
alfa_sigma = 2  #outside good sd

# Three/two products, 3 vars (2 continuous,1 dummy)
threeprod <- list()
twoprod <- list()
purchase <- list()

for (t in 1:1000){
  threeprod[[t]] = cbind(rep(t,3),c(1,1,1),c(1,2,3),runif(3),runif(3),ceiling(runif(3,-0.5,0.5)))
  purchase[[t]] = which.max(rbind(threeprod[[t]][,c(4,5,6)]%*%mvrnorm(1,beta_mu,beta_sigma) + 
    matrix( -log(-log(runif(3))), 3, 1),rnorm(1,alfa_mu,alfa_sigma)) )
  threeprod[[t]] = cbind(threeprod[[t]],c(purchase[[t]],0,0))
}

for (t in 1001:2000){
  twoprod[[t]] = cbind(rep(t,2),c(1,1),c(1,2),runif(2),runif(2),ceiling(runif(2,-0.5,0.5)))
  purchase[[t]] = which.max(rbind(twoprod[[t]][,c(4,5,6)]%*%mvrnorm(1,beta_mu,beta_sigma) + 
    matrix( -log(-log(runif(2))), 2, 1),rnorm(1,alfa_mu,alfa_sigma)) )
  if (purchase[[t]] == 3) {purchase[[t]] <- 4}
  twoprod[[t]] = cbind(twoprod[[t]],c(purchase[[t]],0))
}

X <- rbind(do.call(rbind,threeprod),do.call(rbind,twoprod))

xcoding <- c(1,1,1)

mcmc = list(R = 5000, use = 2000)
options = list(none=TRUE, save=TRUE, keep=5)

out = choicemodelr(X, xcoding, mcmc = mcmc,options = options)