Inverse Predictive Posterior Suppose I have a parametric nonlinear model, say
$$
y_i |\theta \sim N(f_{\theta}(x_i), \sigma^2)
$$
with known form of $f_\theta$. We get data $d=(y_i,x_i)_{i=1,\ldots,n}$ and obtain posterior samples so we can make inference on $\theta|d$. This data is collected in such a way that $x$ is not random, but rather fixed by experimental design.
Now suppose I get new data $y^{(new)}$ but no $x^{(new)}$. Is there someway I can utilize my previous Bayesian model to make posterior inference about $x^{(new)}$?
EDIT: apparently in frequentist world this is called inverse prediction and the initial dataset is called a calibration dataset.
EDIT2: What would the consequences be if I were to fit the inverse model
$$
x_i|\theta \sim N(f^{-1}_\theta(y_i),\sigma_x^2)
$$
Note that $x_i$ in the calibration set are fixed and measured without error, so this seems strange to me and probably problematic in some way. I then could get my answer by using the predictive posterior $p(x^{(new)}|x)$. In retrospect this approach seems nonsensical.
EDIT3: My particular model can also be formulated as $y_i = f_\theta(x_i) + \epsilon_i$ where $\epsilon_i \sim N(0,\sigma^2)$. I could use my posterior samples to draw from this residual $e_j \sim \epsilon|d$ and then $f^{-1}_\theta(y^{(new)} - e_j)$ might be samples from the distribution I am interested in. Thoughts?

Ok, I've edited my response taking into account feedback from the OP.  Below is a DAG that captures the assumptions provided. So for example, $x^{(\mathrm{new})}$ need not be equal in distribution to $x$ as required, and $y^{(\mathrm{new})}$ is conditionally independent of the training data $\mathbf{Y}$ given $\theta$

We require draws from the distribution $p(x^{(\mathrm{new})}|y^{(\mathrm{new})}, \mathbf{Y})$. We have that
\begin{align}
p(x^{(\mathrm{new})}|y^{(\mathrm{new})}, \mathbf{Y}) \propto p(y^{(\mathrm{new})}|x^{(\mathrm{new})}, \mathbf{Y}) p(x^{(\mathrm{new})}|\mathbf{Y})
\end{align}
where we've dropped terms that do not depend on $x^{(\mathrm{new})}$. Focusing on one term at a time:
\begin{align}
p(y^{(\mathrm{new})}|x^{(\mathrm{new})}, \mathbf{Y}) &= \int_\theta p(y^{(\mathrm{new})}|x^{(\mathrm{new})}, \mathbf{Y}, \theta) p(\theta|x^{(\mathrm{new})}, \mathbf{Y})d\theta\\
&= \int_\theta p(y^{(\mathrm{new})}|x^{(\mathrm{new})},\theta) p(\theta|\mathbf{Y})d\theta\\
&= E_{\theta|\mathbf{Y}} \left[p(y^{(\mathrm{new})}|x^{(\mathrm{new})},\theta) \right]
\end{align}
The implied conditional independences from the DAG justify dropping the terms in the second line. Next
\begin{align}
p(x^{(\mathrm{new})}|\mathbf{Y}) \propto p(x^{(\mathrm{new})},\mathbf{Y})\propto p(x^{(\mathrm{new})})
\end{align}
again using the DAG as justification that $x^{(\mathrm{new})}$ and $\mathbf{Y}$ are independent.
Thus,
\begin{align}
p(x^{(\mathrm{new})}|y^{(\mathrm{new})}, \mathbf{Y}) \propto p(x^{(\mathrm{new})}) E_{\theta|\mathbf{Y}} \left[p(y^{(\mathrm{new})}|x^{(\mathrm{new})},\theta) \right]
\end{align}
Then use an importance sampling approach to draw from this distribution. Specifically, first create $M$ draws from the prior distribution for $x$, as in $x^{(\mathrm{new})}_\ell \sim p(x^{(\mathrm{new})})$, $\ell=1,\ldots,M$, where $M$ is some very large integer.
Then calculate $w_\ell\equiv E_{\theta|\mathbf{Y}} \left[p(y^{(\mathrm{new})}|x^{(\mathrm{new})}_\ell,\theta) \right]$ and resample with replacement from your set of $M$ values, with sampling probability proportional to $w_\ell$. I do not know how large your resampled set should be, but I think it should be smaller than $M$.
Your resampled set can be taken as a draw from the posterior predictive distribution of interest.
EDIT: Or do a Metropolis-Hasting algorithm to do your sampling. Or others. Xi'an (who commented on your question) is an expert on sampling algorithms.
