In this paper, we propose an approach to modeling functional magnetic resonance imaging (fMRI) data that combines hierarchical polynomial models, Bayes estimation, and clustering. A cubic polynomial is used to fit the voxel time courses of event-related design experiments. The coefficients of the polynomials are estimated by Bayes estimation, in a two-level hierarchical model, which allows us to borrow strength from all voxels. The voxel-specific Bayes polynomial coefficients are then transformed to the times and magnitudes of the minimum and maximum points on the hemodynamic response curve, which are in turn used to classify the voxels as being activated or not. The procedure is demonstrated on real data from an event-related design experiment of visually guided saccades and shown to be an effective alternative to existing methods.