I can provide an analytical framework for the problem, but I'm not familiar enough with PyMC3 or Stan to translate that framework into either of those languages. (I write my own MCMC samplers in Mathematica!)
Let $y = (y_1, \ldots, y_n)$ denote an unobserved variable where
\begin{equation}
p(y|\mu,\tau) = \prod_{i=1}^n \textsf{N}(y_i|\mu,\tau)
\end{equation}
and $\tau$ is a variance.
The actual observations are given by $z = (z_1, \ldots, z_n)$, where
\begin{equation}
p(z|y,\lambda_1,\lambda_2) = \prod_{i=1}^{n_1} \textsf{N}(z_i|y_i,\lambda_1) \prod_{i=n_1+1}^n \textsf{N}(z_i|y_i,\lambda_2)
\end{equation}
and $\lambda_j$ is a variance. The model is completed by providing a prior for the parameters $\phi = (\mu,\tau,\lambda_1,\lambda_2)$. This prior will likely impose some independence assumptions such as
\begin{equation}
p(\phi) = p(\mu,\tau)\,p(\lambda_1)\,p(\lambda_2) .
\end{equation}
The joint posterior for the unknowns given the observations $z$ can be expressed as
\begin{equation}
p(y,\phi|z) \propto p(z|y,\lambda_1,\lambda_2)\,p(y|\mu,\tau)\,p(\phi) .
\end{equation}
The goal might be to compute the posterior predictive distribution for the latent variable:
\begin{equation}
p(y_{n+1}|z) = \iint \textsf{N}(y_{n+1}|\mu,\tau)\,p(\mu,\tau|z)\,d\mu\,d\tau ,
\end{equation}
where
\begin{equation}
p(\mu,\tau|z) = \iiint p(y,\phi|z)\,dy\,d\lambda_1\,d\lambda_2 .
\end{equation}
Let $\{(y^{(r)},\phi^{(r)}\}_{r=1}^R$ denote draws from the joint posterior distribution (computed via MCMC). The predictive distribution is approximated by
\begin{equation}
p(y_{n+1}|z) \approx \frac{1}{R} \sum_{r=1}^R \textsf{N}(y_{n+1}|\mu^{(r)},\tau^{(r)}) .
\end{equation}
Whether this predictive distribution is approximately normal or not depends strongly on the prior for $\mu$ and $\tau$. For example, if the prior is parametric, such as
\begin{equation}
p(\mu,\tau) = p(\mu)\,p(\tau) = \textsf{N}(\mu|m_0,s_0^2)\,\textsf{Inv-Gamma}(\tau|a_y/2,b_y/2) ,
\end{equation}
then the predictive distribution will likely be approximately normal. By contrast, if one adopts a nonparametric prior for $(\mu,\tau)$ such as the Dirichlet Process (DP) prior, then the posterior predictive distribution will reflect the underlying distribution for $y$, which could be highly non-normal.
The OP asked for more detail on the DP-related prior. As a preliminary, it is convenient to re-express the prior for $y$ as follows:
\begin{align}\label{like_mix}
y_i|\theta_y &\stackrel{\text{iid}}{\sim} F(\theta_y) \\
\label{parametric_H}
\theta_y &\sim H ,
\end{align}
where $\theta_y = (\mu,\tau)$ and the densities for the distributions $F(\theta_y)$ and $H$ are given by
\begin{align}
f(y_i|\theta_y) &= \textsf{N}(y_i|\mu,\tau) \\
h(\theta_y) &= \textsf{N}(\mu|m_0,s_0^2)\, \textsf{Inv-Gamma}(\tau|a_y/2,b_y/2) .
\end{align}
As we acquire data, we learn about $\theta_y$. We now turn to a more general prior that has the property that --- instead of learning about $\theta_y$ --- we learn about the distribution for $\theta_y$.
We generalize the parametric prior by providing a hierarchical prior for $\theta_y$ that involves a Dirichlet Process (DP):
\begin{align}
\theta_y|G &\sim G \\
G &\sim \textsf{DP}(\xi,H) ,
\end{align}
where $\textsf{DP}(\xi,H)$ represents a Dirichlet Process.
In this prior, the prior for $\theta_y$ is the probability measure $G$ which itself is random and depends on $\xi$ and $H$.
The mean of $G$ is $H$ (the base distribution): $E[G] = H$ and the variation of $G$ around $H$ is controlled by the $\xi$ (the concentration parameter).
There is an explicit representation for the DP prior that involves an infinite-order mixture. In preparation, let $\psi = (w,\theta)$ denote the parameters of the mixture where $w = (w_1, w_2, \ldots )$ is an infinite collection of nonnegative mixture weights that sum to one, $\theta = (\theta_1, \theta_2, \ldots)$ is a corresponding collection of mixture parameters, and $\theta_c = (\mu_c,\tau_c)$.
In addition, let $\delta_{\theta_c}$ denote a point mass located at $\theta_c$.
The prior can be expressed as
\begin{align}
\theta_y|\psi &\sim G = \sum_{c=1}^\infty w_c\,\delta_{\theta_c} \\
w|\xi &\sim \textsf{Stick}(\xi) \\
\theta_c &\stackrel{\text{iid}}{\sim} H .
\end{align}
In this representation, the randomness of $G$ follows from the randomness of the mixture weights $w$ (which depend on $\xi$) and of the locations $\theta$ (which depend on $H$). As we acquire data, we learn about $\psi = (w,\theta)$ and consequently we learn about the distribution for $\theta_y$.
The form of the prior that will prove most useful is obtained by integrating out $\theta_y$, which has the effect of transfering the infinite mixture from $\theta_y$ to $y$:
\begin{align}\label{x_mix}
y_i|\psi &\stackrel{\text{iid}}{\sim} \sum_{c=1}^\infty w_c\,F(\theta_c) \\
w|\xi &\sim \textsf{Stick}(\xi) \\
\theta_c &\stackrel{\text{iid}}{\sim} H .
\end{align}
All that remains is to describe $\textsf{Stick}(\xi)$, the distribution for the mixture weights. (In addition, the DPM can be extended by providing a prior for $\xi$ as well.)
The standard prior for $w$ is given by
\begin{equation}
w|\xi \sim \textsf{Stick}(\xi),
\end{equation}
where $\textsf{Stick}(\xi)$ denotes the stick-breaking distribution given by
\begin{equation}\label{gem1}
w_c = v_c \prod_{\ell=1}^{c-1} (1-v_\ell)
\qquad\text{where $v_c \stackrel{\text{iid}}{\sim} \textsf{Beta}(1, \xi)$}.
\end{equation}
The concentration parameter $\xi$ controls the rate at which the weights decline on average. In particular, the weights decline geometrically in expectation:
\begin{equation}
E[w_c|\xi] = \xi^{c-1}\, (1+\xi)^{-c}.
\end{equation}
Note $E[w_1|\xi] = 1/(1+\xi)$ and $E\!\left[\sum_{c=m+1}^\infty w_c|\xi\right] = \big(\xi/(1+\xi)\big)^m$.
