Применение nlminb к подмножествам данных (по индексу или метке) и сохранение того, что программа возвращает как новый фрейм данных - PullRequest
1 голос
/ 02 марта 2010

Мне было интересно, сможет ли кто-нибудь помочь мне в этой, казалось бы, легкой задаче. Я использую nlminb для проведения оптимизации и вычисления статистики по индексу. Вот пример из справки nlminb.

> x <- rnbinom(100, mu = 10, size = 10)
> hdev <- function(par) {
+     -sum(dnbinom(x, mu = par[1], size = par[2], log = TRUE))
+ }
> nlminb(c(9, 12), hdev)
$par
[1] 9.730000 5.954936
$objective
[1] 297.2074
$convergence
[1] 0
$message
[1] "relative convergence (4)"
$iterations
[1] 10
$evaluations
function gradient
      12       27

Предположим, я генерирую случайные величины x, y и z, где z действует как индекс (от 1 до 3).

> x <- rnbinom(100, mu = 10, size = 10)
> y <- rnbinom(100, mu = 10, size = 10)
> z <- rep(1:3, length=100)
> A <- cbind(x,y,z)
> hdev <- function(par) {
+     -sum(dnbinom(x+y, mu = par[1], size = par[2], log = TRUE))}

1) Как я могу применить nlminb(c(9, 12), hdev) к данным, установленным индексом z? Другими словами, я бы хотел вычислить nlminb(c(9, 12), hdev) для z=1, z=2 и z=3 отдельно. Я пробовал by(A, z, function(A) nlminb(c(9,12), hdev)) и sparseby(A, z, function(A) nlminb(c(9,12), hdev)), но они возвращают абсолютно одинаковые значения для каждого значения z.

2) Я хотел бы превратить каждый вывод в новый фрейм данных, чтобы он стал матрицей 3X2.

[1] Z1_ANSWER_1 Z1_ANSWER_2
[2] Z2_ANSWER_1 Z2_ANSWER_2
[3] Z3_ANSWER_1 Z3_ANSWER_2

Поскольку nlminb возвращает сводку статистики, мне нужно было использовать CASEZ1<-nlminb$par, CASEZ2<-nlminb$par, CASEZ3<-nlminb$par, а затем использовать cbind для их объединения. Однако я хотел бы автоматизировать этот процесс, поскольку реальные данные, над которыми я работаю, имеют гораздо больше категорий, чем z, представленное здесь.

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

Заранее большое спасибо.

1 Ответ

1 голос
/ 02 марта 2010

Позвольте мне попробовать подход

x <- rnbinom(100, mu = 10, size = 10)
y <- rnbinom(100, mu = 10, size = 10)
z <- rep(1:3, length=100)
A <- as.data.frame(cbind(x,y,z))

При первой загрузке plyr библиотека

library(plyr)

Следующий код возвращает результаты для каждого z

dlply(A, .(z), function(x) {
    hdev <- function(par, mydata) {-sum(dnbinom(mydata, mu = par[1], size = par[2], log = TRUE))}
    nlminb(c(9, 12), hdev, mydata=t(as.vector(x[1] + as.vector(x[2]))))
}
)

Теперь с этим вы получите фрейм данных 3x2 с результатами $ par

ddply(A, .(z), function(x) {
    hdev <- function(par, mydata) {-sum(dnbinom(mydata, mu = par[1], size = par[2], log = TRUE))}
    res <- nlminb(c(9, 12), hdev, mydata=t(as.vector(x[1] + as.vector(x[2]))))
    return(res$par)
}
)
...