We have used an importance-sampling Monte Carlo technique to study the Baxter-Wu model with various fractions of the lattice sites occupied by random, quenched, nonmagnetic, site impurities. We found the system had long-time fluctuations which were caused by the creation and motion of domain boundaries. For this reason our study required the use of a large number of Monte Carlo steps per spin. The data were analyzed using finite-size scaling to extract the infinite-lattice critical exponents and critical amplitudes from the Monte Carlo simulation of finite lattices. This analysis of the data for the pure Baxter-Wu model yielded results which agreed well with exact results and series-expansion predictions. The addition of quenched site impurities caused a dramatic change in the critical behavior. The pure-lattice critical exponents ($\ensuremath{\nu}=\frac{2}{3}$, $\ensuremath{\alpha}=\frac{2}{3}$, ${\ensuremath{\gamma}}_{m}\ensuremath{\simeq}1.17$) change upon the addition of only a few percent impurities to $\ensuremath{\nu}=1.00\ifmmode\pm\else\textpm\fi{}0.07$, $\ensuremath{\alpha}\ensuremath{\lesssim}0.0$, and ${\ensuremath{\gamma}}_{m}=1.95\ifmmode\pm\else\textpm\fi{}0.08$. The phase diagram as a function of impurity concentration and an estimate for the infinite-lattice percolation limit are also given.