min_rep_fg.R

## You can show that the estimate of fg is wrong *only*
## when m_beta = 1 and m_c = 1 and xx = 1 and that in
## such a case it appears to overcount by 2, implying term 11 in the derivation is wrong

m_alpha <- 0
s_alpha <- 1
m_beta <- 1
s_beta <- 1
m_t <- 0
s_t <- 1
y <- 0
m_c <- 1
s_c <- 1
xx <- 1
alpha_sum <- m_alpha * xx; beta_sum <- m_beta * xx
alpha_square_sum <- xx^2 * (m_alpha^2 + s_alpha)
beta_square_sum <- xx^2 * (m_beta^2 + s_beta)
mu2 <- function() {
  rnorm(1, m_alpha, sqrt(s_alpha)) * xx + 
    rnorm(1, m_c, sqrt(s_c)) * rnorm(1, m_t, sqrt(s_t)) + 
    rnorm(1, m_beta, sqrt(s_beta)) * xx * rnorm(1, m_t, sqrt(s_t))
}
ff <- function() (y - mu2())^2


# g, y, m_t, s_t, m_c, s_c, m_alpha, s_alpha, 
# m_beta, s_beta, m_mu, s_mu, alpha_sum, beta_sum, 
# alpha_square_sum, beta_square_sum)
clvm:::calculate_fg(0, matrix(y), m_t, s_t, m_c, s_c, 
             matrix(m_alpha), matrix(s_alpha), 
             matrix(m_beta), matrix(s_beta), 0, 0, 
             matrix(alpha_sum), matrix(beta_sum), 
             matrix(alpha_square_sum), matrix(beta_square_sum))

fs <- replicate(1e5, ff())
mean(fs)