g_i <- function(x, i, mu, sigma){
  res <- rep(NA, length(x))
  for(k in 1:length(x)){
    res[k] <- x[k]*prod(1-pnorm(x[k], mu[-i], sigma[-i]))*dnorm(x[k],mu[i],sigma[i])
  }
  return(res)
}

expected_score <- function(mu, sigma){
  n <- length(mu)
  res <- sum(mu)
  for(i in 1:n){
    I_i <- integrate(g_i, lower=-Inf, upper=Inf, i=i, mu=mu, sigma=sigma)$value
    res <- res - I_i/sigma[i]
    }
    return(res)
  }