Training a Bernoulli model using probabilities as inputs I'm using two methods to train a Bernoulli model, and am trying to understand why they are not yielding similar results. For both methods, I have a length $N$ array of probabilities $\{\hat{y}^{(n)}\}_{n=1}^{N}$, and I want to estimate the distribution of a length $N$ array of parameters $\{\theta^{(n)}\}_{n=1}^{N}$. In method (1), for each $\hat{y}^{(n)}$, I sample $M$ times from a Bernoulli distribution with probability $\hat{y}^{(n)}$, and use the resulting binary data as the input to my model. In method (2), I use $\hat{y}^{(n)}$ directly by incrementing the log joint density using the following update rule:
$
 \begin{equation*}
  \begin{split}
  \log p(\mathbf{y} \mid \boldsymbol{\theta}) &= \sum_{n=1}^{N} \log p(y^{(n)} \mid \theta^{(n)})\\
  &\approx \sum_{n=1}^{N} \frac{1}{M} \sum_{m=1}^{M} \log p(y^{(n)}_m \mid \theta^{(n)})\\
  &\approx \sum_{n=1}^{N} E_{y^{(n)}_m}[\log p(y^{(n)}_m \mid \theta^{(n)})]\\
  &= \sum_{n=1}^{N} E_{y^{(n)}_m}\left[\log({\theta^{(n)}}^{y^{(n)}_m} \cdot (1 - \theta^{(n)})^{1 - y^{(n)}_m})\right]\\
  &= \sum_{n=1}^{N} E_{y^{(n)}_m}\left[y^{(n)}_m \log(\theta^{(n)}) +  (1 - y^{(n)}_m) \log(1 - \theta^{(n)})\right]\\
  &= \sum_{n=1}^{N} \Pr(y^{(n)}_m = 1) \log(\theta^{(n)}) + \Pr(y^{(n)}_m = 0) \log(1 - \theta^{(n)})\\
  &= \sum_{n=1}^{N} \hat{y}^{(n)} \log(\theta^{(n)}) + (1 - \hat{y}^{(n)}) \log(1 - \theta^{(n)})
  \end{split}
 \end{equation*}
$
I'm using Stan, where you specify the log joint density incrementally using each data point. Pseudocode for these two methods looks like this:

I expect methods (1) and (2) to yield similar estimates for $\theta$ for large $M$, but I'm finding this to not be the case. I've reproduced this issue on a small toy problem using Stan, here is the code:
import matplotlib.pyplot as plt
import numpy as np
import pystan

def get_theta_mean(fit):
    samples = fit.extract()
    theta = np.moveaxis(samples['theta'], 0, -1)
    return theta.mean(axis=1)

rng = np.random.RandomState(0)
N = 100
probs = rng.uniform(0, 1, N)

binary_model = '''
data {
    int<lower=0> N;
    int<lower=0> M;
    int<lower=0, upper=1> y[M, N];
}
parameters {
    real<lower=0, upper=1> theta[N];
}
model {
    for (m in 1:M) {
        y[m] ~ bernoulli(theta);
    }
}
'''
binary_sm = pystan.StanModel(model_code=binary_model)

M_list = [10, 100, 1000]
theta_means = {}

for M in M_list:
    y = np.full((M, N), np.nan)
    for m in range(M):
        for n in range(N):
            y[m, n] = rng.binomial(1, probs[n])
    y = y.astype(int)

    binary_fit = binary_sm.sampling(
        data={'N': N, 
              'M': M, 
              'y': y})

    theta_means[M] = get_theta_mean(binary_fit)

prob_model = '''
data {
    int<lower=0> N;
    real<lower=0, upper=1> yhat[N];
}
parameters {
    real<lower=0, upper=1> theta[N];
}
model {
    for (n in 1:N) {
        target += lmultiply(yhat[n], theta[n]) + lmultiply(1 - yhat[n], 1 - theta[n]);
    }
}
'''
prob_sm = pystan.StanModel(model_code=prob_model)
prob_fit = prob_sm.sampling(
    data={'N': N,
          'yhat': probs})

prob_theta_mean = get_theta_mean(prob_fit)

for M, theta_mean in theta_means.items():
    plt.scatter(theta_mean, probs, label=f'Method (1), M={M}')
plt.scatter(prob_theta_mean, probs, label='Method (2)')
plt.grid(True)
plt.xlabel(r'$E[\theta]$')
plt.ylabel('Probability')
plt.legend()

Here is a scatterplot of the results I get on the toy problem. (1) is supposed to better approximate (2) as $M$ is increased, eventually converging to $E[\theta^{(n)}] = \hat{y}^{(n)}$. (1) seems to follow this, but (2) is significantly off.

UPDATE:
I realized there is an error in line one of my derivation above, and that method (2) actually represents
$
 \begin{equation*}
   \begin{split}
   \log p(\vec{y} \mid \theta) &= \sum_{n=1}^{N} \log p(y^{(n)} \mid \theta^{(n)})\\
   &\approx \sum_{n=1}^{N} \log\left(\frac{1}{M} \sum_{m=1}^{M}  p(y^{(n)}_m \mid \theta^{(n)})\right)
   \end{split}
 \end{equation*}
$
I changed method (1) to reflect this, and now my results between method (1) and (2) are consistent, and neither of them satisfy $E[\theta \mid D] \approx \hat{y}$.

 A: Update My original answer below applies to your original approximation. With respect to your modified approximation, which is different but also a sensible approximation, the issue of scaling by $M$ or not doesn't matter anymore, because then it just becomes a constant shift in the log-posterior (you are essentially deciding whether to add a constant $-\log (M)$). The results, which are the same between Method 1 and Method 2, don't match your intuition because of a shrinkage effect, which I discuss below.
Original answer
Method (2) is giving the expected results; Method (1) is coded incorrectly because it is not scaled by $M$. If you scale each added increment to the log-posterior by $M$, then you will get similar results for large $M$. In more detail:
Method (2) is using a model based upon the continuous Bernoulli distribution. The posterior means you are finding seem to be correct; the issue is that there is very little data and so the results are largely prior-driven. For example, I ran an adapted version of your code for $N=1$ and $\hat y=0.25$. Here is a histogram estimating the posterior distribution of $\theta$, with a red line indicating the value $\hat y = 0.25$:

The posterior mean of $\theta$ is a shrinkage estimator of the observed data $\hat y = 0.25$ and the prior mean 0.5 (things look 'correct' in that case when you have $\hat y = 0.5$ because the data are exactly equal to the prior mean). But the mode of the distribution is approximately at the place you would expect.
As your mathematical derivation suggests, Method (1) should yield approximately equal inferences for large enough $M$. However, by not scaling by $M$, you are artificially creating a great deal of data, which is having the effect of overwhelming the prior (In a comment, I suggested this shouldn't affect point estimates, but I misspoke). See below for two posterior histograms corresponding to Method (1) using $M = 500$, one which properly scales by $M$ and another which does not $M$. You'll observe that the lack of scaling by $M$ creates an extremely tight, and inappropriately confident, distribution around the observed data.


