Пользовательская функция R не дает правильного ответа с помощью dplyr - PullRequest
Новая
       37

Пользовательская функция R не дает правильного ответа с помощью dplyr

0 голосов
/ 25 апреля 2019

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

f_beta <- function(x, elig, par) {
  return (
    dbeta(x, 
          exp(par[1] + par[2] * log(elig)), 
          exp(par[3] + par[4] * log(elig)), 
          log = TRUE))
}

Когда я пытаюсь применить функцию к данным во фрейме данных, она не возвращает правильные результаты. Поля "правдоподобие" и "правдоподобие2" должны возвращать один и тот же результат, но это не так. 

assignments <- df %>%
  select(-cluster) %>%
  crossing(fits) %>%
  mutate(likelihood = prior * dbeta(enrpct, exp(a0 + b0 * log(elig)), exp(a1 + b1 * log(elig)), log = TRUE), 
         likelihood2 = prior * f_beta(enrpct, elig, c(a0, b0, a1, b1))) %>%
  group_by(id) %>%
  top_n(1, likelihood) %>%
  ungroup()

Полный код ниже. 

library(tidyr)
library(dplyr)

# custom beta distribution parameterization
f_beta <- function(x, elig, par) {
  return (dbeta(x, exp(par[1] + par[2] * log(elig)), exp(par[3] + par[4] * log(elig)), log = TRUE))
}

# log-likelihood
ll <- function(x, elig, par) {
  -sum(f_beta(x, elig, par))
}

# optimizer 
fit_beta <- function(x, elig, init = c(-0.5, 0.2, -1.1, 0.3)) {
  m <- optim(par = init, fn = ll, elig = elig, x = x)
  coef <- m$par
  data_frame(a0 = coef[1], b0 = coef[2], a1 = coef[3], b1 = coef[4], number = length(x))
}


####### generate data
n <- 10000
n2 <- 5000

# mixture 1 parameters
a0 <- -1
b0 <- 0.3
a1 <- -2
b1 <- 1

# mixture 2 parameters
a01 <- -1
b01 <- 0.5
a11 <- -1.5
b11 <- 0.8

# generate data 
df <- data.frame(id = 1:n, elig = sample(1:100, size = n, replace = TRUE) * 100)
df$enrpct <- rbeta(n, exp(a0 + b0 * log(df$elig)), exp(a1 + b1 * log(df$elig)))
df2 <- data.frame(id = (n+1):(n+n2), elig = sample(1:100, size = n2, replace = TRUE) * 100)
df2$enrpct <- rbeta(n2, exp(a01 + b01 * log(df2$elig)), exp(a11 + b11 * log(df2$elig)))
df <- rbind(df, df2)

# assign to clusters randomly like in example
df$cluster <- factor(sample(c("A", "B"), nrow(df), replace = TRUE))

# maximization step of E-M algorithm
fits <- df %>%
  group_by(cluster) %>%
  do(fit_beta(.$enrpct, .$elig)) %>%
  ungroup() %>% 
  mutate(prior = number / sum(number))

# expectation step of E-M algorithm
assignments <- df %>%
  select(-cluster) %>%
  crossing(fits) %>%
  mutate(likelihood = prior * dbeta(enrpct, exp(a0 + b0 * log(elig)), exp(a1 + b1 * log(elig)), log = TRUE), 
         likelihood2 = prior * f_beta(enrpct, elig, c(a0, b0, a1, b1))) %>%
  group_by(id) %>%
  top_n(1, likelihood) %>%
  ungroup()

head(assignments)

1 Ответ

2 голосов
/ 25 апреля 2019

Это потому, что вы используете c(a0, b0, a1, b1) при расчете likelihood2.Но каждая из этих переменных является целым столбцом в кадре данных, поэтому использование c() просто объединит их, и вы в конечном итоге будете использовать неправильные значения.

При этом все будет работать: 

f_beta <- function(x, elig, a0, b0, a1, b1) {
  return (dbeta(x, exp(a0 + b0 * log(elig)), exp(a1 + b1 * log(elig)), log = TRUE))
}

likelihood2 = prior * f_beta(enrpct, elig, a0, b0, a1, b1)
...