Bayesian approach to modelling the prediction interval from a parametric model I think I have already answered my question, but am not completely sure and would like some input. Please feel free to correct my terminology or misconceptions, I learned stats by osmosis rather than formal courses so my knowledge is somewhat patchy and I often seem to use the wrong terminology.
We need to predict the future value of a physical quantity, given a set of measurements of this and other quantities at the current time and present past. We have devised an analytic model with a form motivated by physical understanding, with free parameters that have been tuned by training against a set of past data to optimise the prediction accuracy. I've also run the model through a Markov Chain Monte Carlo (MCMC) analysis, so I have all the marginal distributions for the model parameters etc etc. I also tried using machine learning algorithms, in particular random Forest, instead of the analytic model; however for some values of the predictors the data set is sparse and the results from RF do not seem plausible in many cases. The form of the equation we have encodes prior information (physical understanding) which is probably why it performs better.
Now, the model is to be used operationally to make predictions of the values of physical quantities in the future, based on current and previous measurements of the predictors. For this purpose, it is not enough to simply produce a single prediction; we really need to produce a well calibrated prediction interval. In particular, it is clear from the data and expected from theory that the prediction interval will vary considerably for different values of the predictors.
The first point I would make is that I can't see how I can translate the confidence intervals in the parameters into a prediction interval. Ideally I'd utilise the MCMC chain to produce some kind of distribution of the prediction, but that doesn't seem to be the prediction interval. Please correct me if there is a relatively painless way to utilise this information to do what I'm trying to do. 
I've had a crack at working this through from first principles. I'm probably re-inventing the wheel here.
I have examined the residuals of them model prediction compared with the true values and they are sufficiently close to normally distributed. We can easily evaluate the 'overall' prediction interval by simply determining the sample variance of the residuals. However, if we bin up the residuals against the the relevant predictors, we see a clear trend in the variance as a function of the predictor. In a simple case the binning alone is probably good enough to get an empirical handle on the prediction interval, but we want to work out how to do this in general, so we now assume that we have an 'error model' of the same form as the prediction model, and we then fit the parameters of this model against the 'data' (i.e. the residuals).
To accomplish this I have used Bayes' theorem. For an error model given by
$g(\theta,x)$ 
with $\theta$ being model parameters and $x$ the predictors, we have that
$ P(\theta|R) \propto P(R|\sigma) P(\theta) $
With R the model prediction residual. Assuming a uniform prior for $\theta$ and expanding the posterior term we get
$ P(\theta|R) \propto \prod_{i}^n \frac{1}{g(\theta,x_i)} e^{-\frac{R_i^2}{g(\theta,x_i)^2}}$
Using that I can define the log likelihood and do an optimisation or MCMC analysis on the parameters of the error model and therefore arrive at an empirically calibrated error as a function of the predictors. The error model is a shadow of the prediction model; the same form, with the free parameters optimised against the above likelihood.
Operationally, this approach would predict the mean and sigma of a Gaussian, representing a hopefully accurate model for the range of likely values for the quantity being forecast.
So, does this make sense? Are there great flaws in this approach, or existing tools that can do this more easily/accurately?
 A: It seems that your problem calls for Gaussian model with expectation being your function, that you know apriori from the nature of physical phenomena. Suppose, that physics tells you, that having covariate $X$ with measurements $X = (X_{1},...,X_{n})$ (e.g. temperature, pressure, etc.) the response is a function $E[Y|X]=f(X;\theta)$, where $Y$ is a measurements of response, $\theta$ is model parameters and suppose that $Y \in R$. In addition, suppose that your measurements have Gaussian errors with variance $ \sigma^{2}$. I would really suggest you to employ Bayesian modelling, since posterior distributions fully account for uncertainties in your model (even in small samples!!). 
If you also do not know the values of this covariance matrix, you can put a noninformative prior distribution on it, e.g. Jeffreys prior:  
$$\pi \left ( \sigma  \right )\propto \left (\sigma  \right )^{-1}$$
And noninformative priors for regression parameters could be: $$\pi(\theta)\propto1$$
So, we have a model:
$$Y_{i}\sim N\left ( f\left ( X_{i} \right ),\sigma  \right ), i=1,2,...,n$$
$$\pi \left ( \sigma  \right )\propto \left (\sigma  \right )^{-1}$$
From this model, construct posterior:
$$\pi(\theta, \sigma|Y)\propto L(Y;\theta)\pi(\sigma)\pi(\theta)$$
Then obtain samples from this posterior by MCMC. To obtain prediction at specific covariate value $X^{*}$ use this MCMC sample of parameters $\theta$ and $\sigma$ to generate random variables from your Gaussian model. When this sample is obtained you could use its average as estimate of the prediction and the prediction interval can be constructed by constructing empirical cdf. I would really suggest you to use "Bayesian analysis" by Andrew Gelman.
A: One option for finding the prediction interval or the posterior predictive distribution is to use one of the *bugs software packages (winbugs, openbugs, jags, etc.) or something similar, then include an additional datapoint with the value missing.  You will probably need to include a diffuse prior for the missing value.  Then each iteration of the McMC algorithm will estimate the missing value from the likelihood using the current estimates of the parameters.  This takes into account the uncertainty in the parameters and the variability within the distribution of the individual points.  The realizations for the missing value can be used to estimate the posterior predictive distribution for a new observation and an interval (HPD or other) from that distribution will give the prediction interval. 
