When does the marginal MLE converge to the complete data MLE? What I mean by the title is suppose we have a distribution $p(x,z\;|\;\theta)$, where the $x$ are observed and each $x_i$ depends on a hidden $z_i$.  Then the marginal MLE is given by
$$\max_\theta\prod_{i=1}^N\int p(x_i,z\;|\;\theta)dz$$
If we take $N\rightarrow\infty$, under what circumstances will this converge to the true $\theta$?  That is to say, under what circumstances will it converge to the $\theta$ the complete MLE would converge to when the complete MLE is derived from the complete-data case where the $z_i$ are known?
The example I'm thinking of is the coin-flipping example, where we choose coin $A$ or coin $B$ with equal probability, and then flip the coin $F$ times.  So the marginal likelihood would be
$$\max_{\theta}\prod_{i=1}^N\frac{1}{2}\binom{F}{x_i}\big[\theta_A^{x_i}(1-\theta_A)^{F-x_i}+\theta_B^{x_i}(1-\theta_B)^{F-x_i}\big]$$
I'm quite confident that if $F=1$ then even as $N\rightarrow\infty$ our estimate for $\theta$ will never converge to the true $\theta$, but if $F>1$ I'm thinking maybe it will, but I'm not sure.
 A: A very pertinent question! The general and unsurprising answer is that the marginal will have the properties associated with the density$$p(x|\theta)=\int p(x_i,z\;|\;\theta)\text{d}z$$rather than with the complete density $p(x,z|\theta)$. If $\theta$ is identifiable in the marginal model, i.e., for the density $p(\cdot|\theta)$, then the MLE should converge to the true $\theta$ as $n$ grows to infinity, under the usual provisions for the regularity of the density, &tc.
However, the MLE will not have the same convergence properties as the complete data MLE in that, e.g., its asymptotic variance will not be the same, given that the Fisher information is smaller for the marginal. See for instance the reference book by Rubin and Little.
In the mixture models, the parameter may be identifiable up to a permutation of the labels, in which case MLEs may converge when well-defined. For instance, in a Gaussian mixture with all parameters unknown the likelihood is unbounded and hence the MLE is not defined stricto sensu. However, there exists a local mode converging to the true value of the parameter at the usual speed of $\sqrt{N}$. See, e.g., the book by McLachlan and Peel.

In your example of a mixture of two Bernoulli distributions, the parameter $(\theta_A,\theta_B)$ is not identifiable, given that the mixture of two Bernoullis is again a Bernoulli with parameter $\{\theta_A+\theta_B\}/2$. As shown by the picture above, obtained from 10⁵ observations from a mixture with $(\theta_A,\theta_B)=(.1,.5)$, the marginal log-likelihood is constant along the lines $\{\theta_A+\theta_B\}/2=c$. While the complete log-likelihood identifies the true $(\theta_A,\theta_B)$.

Now, if you move to $F=2$, i.e., to a mixture of two Binomial $\mathfrak{B}(F,\theta_A)$ and $\mathfrak{B}(F,\theta_B)$, the parameter $(\theta_A,\theta_B)$ is identifiable, up to a permutation of the probabilities, as shown by the picture above, associated with 5555 observations from a Binomial mixture with $(\theta_A,\theta_B)=(.1,.5)$
P.S.: The graphs were obtained with the following R code:
theta=c(.1,.5)
n=1e5
compo=sample(c(1,2),n,rep=TRUE)
dat=as.integer((runif(n)<theta[compo]))
mdat=1-dat

marlik=function(theta){
   sum(log(theta[1]^dat*(1-theta[1])^mdat+theta[2]^dat*(1-theta[2])^mdat))}

complik=function(theta){
   sum(log(theta[compo])*dat+log(1-theta[compo])*mdat)}

theta1=seq(.01,.99,le=123)
theta2=seq(.01,.99,le=132)

complix=malix=matrix(NA,123,134)
for (i in 1:123)
for (j in 1:132){
malix[i,j]=marlik(c(theta1[i],theta2[j]))
complix[i,j]=complik(c(theta1[i],theta2[j]))}

par(mfrow=c(1,2))
image(malix,col=terrain.colors(99),xlab=expression(theta[1]),ylab=expression(theta[2]));contour(malix,add=TRUE,nlevels=99)
image(complix,col=terrain.colors(99),xlab=expression(theta[1]),ylab=expression(theta[2]));contour(complix,add=TRUE,nlevels=99)

n=5555
dat=as.integer((runif(n)<theta[compo]))+as.integer((runif(n)<theta[compo]))
compo=sample(c(1,2),n,rep=TRUE)

marlik=function(theta){
   sum(log(dbinom(dat,2,theta[1])+dbinom(dat,2,theta[2])))}

complik=function(theta){
   sum(dbinom(dat,2,theta[compo],log=TRUE))}

complix=malix=matrix(NA,123,134)
for (i in 1:123)
for (j in 1:132){
malix[i,j]=marlik(c(theta1[i],theta2[j]))
complix[i,j]=complik(c(theta1[i],theta2[j]))}

par(mfrow=c(1,2))
image(malix,col=terrain.colors(99),xlab=expression(theta[1]),ylab=expression(theta[2]));contour(malix,add=TRUE,nlevels=99)
image(complix,col=terrain.colors(99),xlab=expression(theta[1]),ylab=expression(theta[2]));contour(complix,add=TRUE,nlevels=99)

