4

トゥイーディー (または複合ポアソンガンマ) 分布の有限混合を推定しようとしています。これを行う方法に関するリソースを見つけることなく、考えられるリソースを精査しました。

現在、R でflexmixパッケージを使用して、 12 ～ 14 ページのflexmixビネットで概説されているように、別の M ステップ ドライバーを作成しようとしています。cplmパッケージに依存する私のコードは次のとおりです。

tweedieClust <- function(formula = .~.,offset = NULL){
require(tweedie)
require(cplm)
require(plyr)
require(dplyr)
retval <- new("FLXMC", weighted = TRUE, formula = formula, dist = "tweedie",
              name = "Compound Poisson Clustering")

retval@defineComponent <- expression ({
    predict <- function(x, ...) {
        pr <- mu
    }
    logLik <- function(x, y, ...){
        dtweedie(y, xi = p, mu = mu, phi = phi) %>%
             log
    }
    new("FLXcomponent",
        parameters=list(coef=coef),
        logLik=logLik, predict=predict,
        df=df)
})
retval@fit <- function (x, y, w, component) {
    fit <- cpglm(formula = y ~ x, link = "log", weights=w, offset=offset)
    with(list(coef = coef(fit), df = ncol(x),mu = fit$fitted.values,
              p = fit$p, phi = fit$phi),
         eval(retval@defineComponent))
}
retval
}

ただし、これにより次のエラーが発生します。

dtweedie(y, xi = p, mu = mu, phi = phi) のエラー: 非適合配列でのバイナリ演算

tweedie ディストリビューションの有限混合を実行または見た人はいますか? flexmixなどを使用して、これを達成するための正しい方向を教えてもらえますか?

4

1 に答える 1

2

問題は重み部分のどこかにあり、それを削除すると機能します。

tweedieClust <- function(formula = .~.,offset = NULL){
  require(tweedie)
  require(statmod)
  require(cplm)
  require(plyr)
  require(dplyr)
  retval <- new("FLXMC", weighted = F, formula = formula, dist = "tweedie",
            name = "Compound Poisson Clustering")

  retval@defineComponent <- expression ({
    predict <- function(x, ...) {
      pr <- mu
    }
    logLik <- function(x, y, ...){
      dtweedie(y, xi = p, mu = mu, phi = phi) %>%
        log
    }
    new("FLXcomponent",
        parameters=list(mu=mu,xi=p,phi=phi),
        logLik=logLik, predict=predict,
        df=df)
  })
  retval@fit <- function (x, y, w, component) {
    fit <- cpglm(formula = End~.,data=dmft, link = "log")
    with(list(df = ncol(x), mu = fit$fitted.values,
              p = fit$p, phi = fit$phi),
         eval(retval@defineComponent))
  }
  retval
}

例：

library(flexmix)
data("dmft", package = "flexmix")
m1 <- flexmix(End ~ .,data=dmft, k = 4, model = tweedieClust())
于 2015-10-12T07:50:55.170 に答える