EM algorithm for normal mixtures with constraints I have $G$ groups, each with $N_g$ data points $y_{ig}$, $g=1,\dots,G$ and $i = 1, \dots, N_g$. The group for each data point in observable.
I want to estimate the normal mixtures model with $K$ unobservable types, assuming that the type proportions $q_{kg}, k = 1,\dots, K$ are common across groups; i.e. for each $k$ I assume $$q_{k1} = q_{k2} = \dots = q_{kG} = q_k, \quad \text{where} \quad \sum_k q_k = 1. \tag1$$
In other words, I assume a data point $y_{ig}$ is drawn from a normal distribution $\Phi(\cdot; \mu_{kg}, \sigma_{kg})$ with probability $q_{kg} = q_k$. The problem is, how to estimate the parameters $\mu_{kg}, \sigma_{kg}$ and $q_{k}$.
By using expectation-maximization algorithm, I could easily estimate the model separately for each group, but the type proportions $q_{kg}$ would then not agree across different groups. Thus, I need to implement the constraint (1) in the estimation procedure. The question is, how to do it in the context of the EM algorithm?
I use R for the numerics, although as a first step a mathematical explanation suffices.
Edit: an example
Here is an example I have tried to solve: First, I simulate the data, and then I try to estimate it by using the EM algorithm, but something goes wrong, and I am not able to reproduce all the parameters from the simulated data set.
I assume I observe each individual's group $g \in \{1,2\}$ together with three iid data points $y_{it} \sim \Phi(\cdot; \mu_{kg}, \sigma_{kg}), t \in \{1,2,3\},$ for each individual. Here $k \in \{1,2,3\}$ is the unobservable latent type of the individual.
require(data.table)

lognormpdf <- function(Y,mu=0,sigma=1)  {
  -0.5 * (  (Y-mu) / sigma )^2   - 0.5 * log(2.0*pi) - log(sigma)  
}

logsumexp <- function(v) {
  vm = max(v)
  log(sum(exp(v-vm))) + vm
}

logrowSumexp <- function(v) {
  vm = max(v)
  log(rowSums(exp(v-vm))) + vm
}

model.fun <-function(nk) {
  
  model = list()
  # model for Y1,Y2,Y3|k 
  model$A     = array(3*(1 + 3*runif(3*nk)),c(nk))
  model$S     = array(1,c(nk))
  model$pk    = c(0.1,0.5,0.4) # rdirichlet(1,rep(1,nk))
  model$nk    = nk
  return(model)
}


model.simulate <-function(model,N,sd.scale=1) {
  
  Y = array(0,sum(N)) 
  K  = array(0,sum(N)) 
  
  A   = model$A
  S   = model$S
  pk  = model$pk
  nk  = model$nk
  
  K = sample.int(nk,N,TRUE,pk)
  
  Y1  = A[K] + S[K] * rnorm(N) *sd.scale
  Y2  = A[K] + S[K] * rnorm(N) *sd.scale
  Y3  = A[K] + S[K] * rnorm(N) *sd.scale
  
  data.sim = data.table(y1=Y1,y2=Y2,y3=Y3)
  
  return(data.sim)  
}

N = 100000
nk = 3
ng = 2

set.seed(123432)

model = model.fun(nk)
model2 = model.fun(nk)
data = model.simulate(model, N/2)
data2 = model.simulate(model2, N/2)

data$g = 1
data2$g = 2

data = rbind(data, data2)

Y1 = data$y1
Y2 = data$y2
Y3 = data$y3
G = data$g


tau = array(0,c(N,nk))
lpm = array(0,c(N,nk))
lik0 = 1
lik = 0
pk = array(c(1/3,1/3,1/3),c(3,1))



A = array(rnorm(9, 5, 1), c(nk,ng))
S = array(abs(rnorm(9, 1, 10)), c(nk,ng))
tol = 1e-24
iter = 0
# tauold = 1


Xg = array(0, c(N, ng))
Xg[,1] = as.numeric(G ==1)
Xg[,2] = as.numeric(G ==2)



while (abs(lik-lik0)>tol) {
  lik0 = lik
  lik = 0
  
  ###### E-step #########
  
  ltau = array(1,N)%*%t(log(pk))
  lnorm1 = lognormpdf(Y1, t(A[,G]), t(S[,G]))
  lnorm2 = lognormpdf(Y2, t(A[,G]), t(S[,G]))
  lnorm3 = lognormpdf(Y3, t(A[,G]), t(S[,G]))
  lall = ltau +lnorm1 + lnorm2 +lnorm3
  lik = logsumexp(lall)
  
  
  # Q1 = sum(tau*lpm)
  # Q2 = sum(tau*lall)
  
  tau = exp(lall - logrowSumexp(lall))
  lpm = lall
  
  # H1 = -sum(tauold*log(tauold))
  # H2 = -sum(tauold*log(tau))
  # tauold = tau

  
  
  pk = colMeans(tau)
  
  
  ##### M-step ########
  
  
  # Solve for means by group
  
  A = t(tau* (Y1 %*% t(rep(1,nk))+Y2 %*% t(rep(1,nk))+Y3 %*% t(rep(1,nk)))) %*% Xg/(3*t(tau) %*% Xg)
  
  # Solve for standard deviations by group
  
  S = sqrt((t(tau* ((Y1 %*%t(rep(1,3))-t(A[,G]))^2+(Y2 %*%t(rep(1,3))-t(A[,G]))^2+(Y3 %*%t(rep(1,3))-t(A[,G]))^2)) %*% Xg)/(3*t(tau) %*% Xg))
  

  # print(c(iter, abs(lik-lik0), lik, Q2-Q1, H2-H1))
  print(c(iter, abs(lik-lik0), lik))
  iter = iter+1
}

 A: If the observed log-likelihood reads
$$\sum_{g=1}^G \sum_{i=1}^{N_g} \log \sum_{k=1}^K q_k \varphi(y_{ig};\mu_{gk},\sigma_{gk})$$
the completed log-likelihood reads
$$\sum_{g=1}^G \sum_{i=1}^{N_g} \sum_{k=1}^K z_{kig}\log \{q_k \varphi(y_{ig};\mu_{gk},\sigma_{gk})\}$$
where the $z_{kig}$'s are the component indicator variables. The E-function can thus be written as
\begin{align}Q(\theta_0,\theta)&=\mathbb E_{\theta_0}\left[\sum_{g=1}^G \sum_{i=1}^{N_g} \sum_{k=1}^K Z_{kig}\log \{q_k \varphi(y_{ig};\mu_{gk},\sigma_{gk})\} \big| \mathbf y\right]\\
&=\sum_{g=1}^G \sum_{i=1}^{N_g} \sum_{k=1}^K \mathbb E_{\theta_0}\left[Z_{kig}\big| \mathbf y\right]\log \{q_k \varphi(y_{ig};\mu_{gk},\sigma_{gk})\}\\
&= \sum_{k=1}^K \log \{q_k\}\sum_{g=1}^G \sum_{i=1}^{N_g}\mathbb P_{\theta_0}\left[Z_{kig}=1\big| \mathbf y\right]\\
&\quad+\sum_{g=1}^G \sum_{i=1}^{N_g} \sum_{k=1}^K \mathbb P_{\theta_0}\left[Z_{kig}=1\big| \mathbf y\right]\log \{\varphi(y_{ig};\mu_{gk},\sigma_{gk})\}\\
\end{align}
which means that the update of the $q_k$'s in the EM algorithm is based on the expected allocations across groups. Nothing changes for the other group-specific parameters.
