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?

This is an observation rather than an answer, and if $f_{\theta}$ has some heinous form and/or your sample size is small may not have any practical relevance for your particular problem.
You know the form of $f_{\theta}$, and that must help! In particular all plausible potential values of $X_{new}|y_{new}$ must be consistent with what you know about the form of $f_{\theta}$.
Now if $f_{\theta}$ belongs to a family of well-behaved functions then that might tell you a great deal about what are the plausible values of $x_{new}$ once you have observed $y_{new}$. For example, if $f_{\theta}$ is continuous in $x$ or better yet (in order of increasing niceness) Lipschitz continuous, or differentiable, or monotonic, or convex, or linear.
In any case it certainly makes sense to look at observed $(x,y)$ pairs in the data, with similar $y$-values to $y_{new}$. As a very simple toy example, suppose that for all observed $y$-values in the data within an interval $I:= y_{new} \pm \epsilon$ it is true that the corresponding observed $x$-values are tightly clustered in a small region $A$ (e.g. a ball or an interval) of the domain of $X$. Clearly in this situation if $f_{\theta}$ is continuous in $x$ we should reason that $P(x_{new}\in A|y_{new}\in I) \gg P(x_{new}\notin A|y_{new}\in I)$.
