Pitfalls of posterior simulation when analysis didn't begin as Bayesian

I've got a situation where I'd like to evaluate a function of a fitted model, and account for the uncertainty in the fitted model. For example, say I want to calculate the minimum of the function $f(.)$, but I can only get $\hat f(.)$, because I don't know the true function. I got $\hat f(.)$ from a linear model, implicitly treating it in the frequentist way -- there is some true $f(.)$ that I can't observe, and my estimates are random variables, with randomness coming from sampling variability. As such, my variance-covariance matrix (and thus my confidence intervals on parameters) are interpreted as giving an estimate of the size and extent of the region where the true, fixed parameters are most likely to lie, given the random sample that I have. Fine.

However, the multivariate normality of the parameter covariance matrix is attractive for trying to solve my problem with trying to estimate the minimum of $f()$. If I were a Bayesian (and maybe I would be if I had time to really learn the methods completely...) I'd view the parameter covariance matrix not as the variability of the estimator, but as the variability of the parameters. I could then take draws from it, evaluate the minimum of $\hat f()$, and get a credible interval for the minimum.

My question is this: what pitfalls would I face if I simulate parameters from a variance-covariance matrix estimated by frequentist methods? If I "change hats", announce "Now think of this problem like a Bayesian", and go forth with posterior simulation? What would a full-Bayesian do differently from the beginning, assuming zero prior information? And would they get a different answer?

Here is a dummy example in R:

Simulate a quadratic function with some noise, and fit a model:

set.seed(1)
x = sort((runif(1000)*5)-2.5)
plot(density(x))
y = x^2+rnorm(1000,sd=3)
m = lm(y~poly(x,2))

Get the value that minimizes the estimated function:

xgrid=seq(from=-2.5,to=2.5,by=.001)
xg = cbind(1,xgrid,xgrid^2)
predmin = xgrid[which.min(xg%*%m\$coef)]

Simulate from the "posterior" and define a function to get the minima of the simulates

library(MASS)
br = mvrnorm(1000,coef(m),vcov(m))
getmin = function(brep){xgrid[which.min(xg%*%brep)]}

Here is the output ("credible" intervals):

1> predmin
 -0.05
1> quantile(apply(br,1,FUN=getmin),prob=.025)
2.5%
-0.104
1> quantile(apply(br,1,FUN=getmin),prob=.975)
97.5%
0.003025

Would a fully Bayesian treatment from the start get a different answer? Is there an explicitly frequentist way of doing the same analysis?

EDIT: Or would it be sufficient to say "supposing my sample is a random draw from a meta-population, these simulates are the parameters that I'd estimate from those many samples from that meta-population"? If so, I've got a frequentist interpretation for what I'd like to do.