2016-09-06 7 views
0

Я пытаюсь соответствовать 2 компонентной модели гаммы смеси к моим данным (остаточные значения, полученные после выполнения обобщенной линейной модели), используя следующую команду (часть коды):Ошибка в gammamixEM ... Попробуйте другое количество компонентов? в R

expr_mix_gamma <- gammamixEM(expr_glm_residuals, lambda = c(0.75,0.25), k = 2, epsilon = 1e-08, maxit = 1000, maxrestarts=20, verb = TRUE) 

код работает для несколько файлов генов (в цикле). он отлично работает для некоторых файлов, тогда как для других он вызывает следующую ошибку:

Error in gammamixEM(expr_glm_residuals, lambda = c(0.75, 0.25), k = 2, : Try different number of components?  

Я не могу понять, что происходит. может ли кто-нибудь пролить свет на то же самое? Благодаря

ответ

0

Мало мнения экспертов:

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

2) От Мартина Maechler, ETH Zurich: Если это случай смесей один мерных гауссиан, также настоятельно рекомендуется смотреть на CRAN пакета «nor1mix» («1» для " один -мерном). Некоторое время теперь, когда небольшой пакет предоставляет альтернативный к EM, а именно прямому ОМП, просто используя Optim(), где вероятность использует несколько смарт параметризации. конечно, как EM, это также зависит от начального значения, , но мой (ограниченный) опыт был t nor1mix :: norMixMLE() работает значительно быстрее и надежнее, чем EM (который I также предоставляет как nor1mix :: norMixEM().

Аппограмма «Начальная стоимость»: на странице справки показано, как использовать kmeans() для «несколько» надежных запусков; в качестве альтернативы, я бы рекомендовал использовать cluster :: pam(), чтобы начать работу там.

Для меня попытка первого подхода улучшилась для некоторых генов, но все еще продолжалась неудача для других. Итак, наконец, я использовал tryCatch, чтобы запустить мой скрипт и потерять информацию для некоторых генов.