Как рассчитать плотность вероятности для отдельных и для комбинированных компонентов из me.weighted ()

Я использую Mclust для оценки вероятности членства в компонентах, но плотность не включается в вывод me.weighted (). Следовательно, я не могу изобразить плотность вероятности. Следующий код является длинным, потому что я хочу четко проиллюстрировать мою цель и проблему, но я четко указываю, где возникает моя проблема / вопрос. Моя последняя часть кода - это моя попытка решения, но она, вероятно, только подчеркивает мое незнание плотностей вероятностей.

В этом исследовательском проекте моя первая цель - вычислить индекс численности рыбы возраста 1 для последующего анализа. Для этого я хочу оценить долю рыбы возраста 1 определенной длины (т.е. ключ длины к возрасту). Разумно предположить, что меньшая мода - это в основном рыба возраста 1, а большая мода - рыба возраста 2+. Мои данные - это длина тела рыбы (длина вилки, см) и численность в процентах от общей численности (т.е. взвешенная одномерная величина). Обратите внимание, что некоторые очень большие длины с небольшими пропорциями были опущены; таким образом, сумма (dat.df $ пропорции) ‹1.

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

Я читал соответствующие статьи (Мерфи; Скрукка и др .; Миньян; R-Bloggers и т. Д.), Но не нашел ответа.

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

Пакеты

library(ggplot2)
library(mclust) 

Данные

dat.df <- data.frame(flcm = 15:33, proportion = c(0.0043, 0.0114, 0.0296, 0.0519, 0.0540, 0.0403, 0.0294, 0.0152, 0.0257, 0.0793, 0.1458, 0.1505, 0.1277, 0.0909, 0.0389, 0.0308, 0.0121, 0.0101, 0.0085), z1 = c(rep(1,9), rep(0,10)), z2 = c(rep(0,9), rep(1,10)))

Данные графика

ggplot()+
  geom_bar(aes(x=dat.df$flcm, y=dat.df$proportion),
           fill = "gray", position="dodge", stat="identity")+
  xlab("Fork length (cm)")+
  ylab("Probability density")+
  theme_bw()

БЕЗ ВЕСА (т.е. игнорировать пропорцию dat.df $)

Подходит для модели смеси без грузов

mod1 <- densityMclust(dat.df[, "flcm"], modelName = "V")

Постройте плотность вероятности

plot(mod1, what = "density", data = dat.df$flcm, breaks = 5)

С ВЕСОМ (т. е. с учетом доли dat.df $)

Установить модель с грузами

mod1_w <- me.weighted(modelName = "V", 
            data = dat.df$flcm,
            z = cbind(dat.df$z1, dat.df$z2),
            weights = dat.df$proportion)

Данные графика с предполагаемым дробным членством (обновлено z)

ggplot()+
  geom_bar(aes(x=dat.df$flcm, y=dat.df$proportion),
           fill = "gray", position="dodge", stat="identity")+
  geom_line(aes(x = dat.df$flcm, 
                y = (mod1_w$z[,1] * dat.df$proportion)), 
            color = "red") +
  geom_line(aes(x = dat.df$flcm, 
                y = (mod1_w$z[,2] * dat.df$proportion)), 
            color = "red") +
    geom_line(aes(x = dat.df$flcm, 
                y = (mod1_w$z[,1] * dat.df$proportion) +
                  mod1_w$z[,2] * dat.df$proportion), 
              color = "green") +
  xlab("Fork length (cm)")+
  ylab("Probability density")+
  theme_bw()

Постройте график плотности вероятности - Вот где возникает моя проблема / вопрос

plot(mod1_w, what = "density", data = dat.df$flcm, breaks = 5)`

Вот моя попытка решения. По сути, для каждого компонента (age1, age2) умножьте вероятности и масштабируйте пропорционально численности:

#age1 probability density
age1 <- mod1_w$z[,1]* #probability of age1 membership multiplied by
  dnorm(dat.df$flcm, mod1_w$parameters$mean[1], #probability of flcm given age1
        mod1_w$parameters$variance$sigmasq[1])* 
  sum(mod1_w$z[,1]*mod1_w$weights) #and scaled to proportional abundance of age1

#age2 probability density
age2 <- mod1_w$z[,2]* #probability of age2 membership multiplied by
  dnorm(dat.df$flcm, mod1_w$parameters$mean[2],
        mod1_w$parameters$variance$sigmasq[2])* #probability of flcm given age2
  sum(mod1_w$z[,2]*mod1_w$weights) #and scaled to proportional abundance of age2

#combined ages probability density
age_all <- age1 + age2

#looks bad - the probability densities don't correspond well with proportional abundance
ggplot()+
  geom_bar(aes(x=dat.df$flcm, y=dat.df$proportion),
           fill = "gray", position="dodge", stat="identity")+
  geom_line(aes(x = dat.df$flcm, 
                y = age1), 
            color = "red") +
  geom_line(aes(x = dat.df$flcm, 
                y = age2), 
            color = "red") +
    geom_line(aes(x = dat.df$flcm, 
                y = age_all), 
              color = "green") +
  xlab("Fork length (cm)")+
  ylab("Probability density")+
  theme_bw()