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?

I'm going to take a Bayesian approach to this.
As far as I can tell, the existence of the training set is irrelevant for this problem -- it doesn't matter how we obtained the model, we can just take it as given and fixed.
So, the actual inference problem is to obtain posterior distribution for $\vec{x}' = [x_1', x_2' \ldots x_N']$ corresponding to an ensemble $\vec{y}' = [y_1', y_2' \ldots y_N']$ (these are the $y^{(new)}$).
Apply Bayes' theorem to the get the conditional distribution
$$
P(\vec{x}' \vert \vec{y}'; \theta) = \frac{P(\vec{y}' \vert \vec{x}'; \theta)P(\vec{x}'; \theta) }{ \int P(\vec{y}' \vert \vec{x}'; \theta) P(\vert \vec{x}'; \theta) }
$$
$$
P(\vec{x}' \vert \vec{y}'; \theta) = \frac{P(\vec{x}'; \theta) \prod_i P(y_i' \vert x_i'; \theta) }{ \int P(\vec{y}' \vert \vec{x}'; \theta) P(\vert \vec{x}'; \theta) }
$$
Now comes the modeling part.  Given the statement "x is not random, but rather fixed by experimental design", we're going to have to make some kind(s) of assumptions about the process that generated the $\vec{y}'$.
The place I'd start is just to assume that the $x_i$ are independently draw from some a priori distribution $P(\vec{x}') = \prod_i P(x_i')$, suitably specified by domain knowledge or by using the relevant maximum entropy prior.  Then the probability distribution fully factorizes, and one can compute the posterior distribution for each factor corresponding to a particular $x_i'$.
If there are domain reasons to assume that the sequence of $x_i$ are correlated with one another, then you'd need to incorporate that into the model.  Note that even in this case, if desired, one can compute $P(x_i' \vert \vec{y}'; \theta)$ by marginalizing out all of the other $x'$ variables.
From here, the main problem is practical: the fact that the factors $P(y_i' \vert x_i'; \theta) = N(y_i' ; f_\theta (x_i') , \sigma^2)$ means that they are not simple functions of $x_i'$ due to the non-linear $f_\theta$ means that you'll probably need to resort to some sort of numerical approach (Monte Carlo maybe), or approximation (saddlepoint method maybe) to construct a representation of the posterior.
