This Ph.D. thesis concerns the application of machine learning methods to hemodynamic models for BOLD fMRI data. Several such models have been proposed by different researchers, and they have in common a basis in physiological knowledge of the hemodynamic processes involved in the generation of the BOLD signal. The BOLD signal is modelled as a non-linear function of underlying, hidden (non-measurable) hemodynamic state variables. The focus of this thesis work has been to develop methods for learning the parameters of such models, both in their traditional formulation, and in a state space formulation. In the latter, noise enters at the level of the hidden states, as well as in the BOLD measurements themselves. A framework has been developed to allow approximate posterior distributions of model parameters to be learned from real fMRI data. This is accomplished with Markov chain Monte Carlo (MCMC) sampling techniques, including 'parallel tempering', an improvement of basic MCMC sampling. On top of this, a method has been developed that allows comparisons to be made of the quality of these models. This is based on prediction of test data, and comparisons of learnt parameters for different training data. This gives estimates of the generalization ability of the models, as well as of their reproducibility. The latter is a measure of the robustness of the learnt parameters to variations in training data. Together, these measures allow informed model comparison, or model choice. Using resampling techniques, a measure of the uncertainty about the generalization ability and reproducibility of the models is also obtained. The results show that for some of the data, the standard so-called 'balloon' model is sufficient. More complex data have also been designed, however, and for these, the stochastic state space version of the standard balloon model is shown to be superior, although an augmented version of the standard balloon model is not found to be an improvement for either data set.