We introduce a new approach to Bayesian pRF model estimation using Markov Chain Monte Carlo (MCMC) sampling for simultaneous estimation of pRF and hemodynamic parameters. To obtain high performance on commonly accessible hardware we present a novel heuristic consisting of interpolation between precomputed responses for predetermined stimuli and a large cross-section of receptive field parameters. We investigate the validity of the proposed approach with respect to MCMC convergence, tuning and biases. We compare different combinations of pRF - Compressive Spatial Summation (CSS), Dumoulin-Wandell (DW) and hemodynamic (5-parameter and 3-parameter Balloon-Windkessel) models within our framework with and without the usage of the new heuristic. We evaluate estimation consistency and log probability across models. We perform as well a comparison of one model with and without lookup table within the RStan framework using its No-U-Turn Sampler. We present accelerated computation of whole-ROI parameters for one subject. Finally, we discuss risks and limitations associated with the usage of the new heuristic as well as the means of resolving them. We found that the new algorithm is a valid sampling approach to joint pRF/hemodynamic parameter estimation and that it exhibits very high performance.