We propose a method for the analysis of functional magnetic-resonance (fMR) data, based on a Bayesian-network representation. Our method identifies multivariate linear/nonlinear voxel-activation pattern differences across groups, which may provide information complementary to that resulting from a general linear model (GLM)-based analysis. In addition, we describe a model-stabilization method based on data resampling, which may be helpful in the presence of small numbers of subjects, or when data are noisy.