  We measure the high-mass stellar initial mass function (IMF) from resolved
stars in M33 young stellar clusters. Leveraging \textit{Hubble Space
Telescope's} high resolving power, we fully model the IMF probabilistically. We
first model the optical CMD of each cluster to constrain its power-law slope
$\Gamma$, marginalized over other cluster parameters in the fit (e.g., cluster
age, mass, and radius). We then probabilistically model the distribution of MF
slopes for a highly strict cluster sample of 9 clusters more massive than
log(Mass/M$_{\odot}$)=3.6; above this mass, all clusters have well-populated
main sequences of massive stars and should have accurate recovery of their MF
slopes, based on extensive tests with artificial clusters. We find the ensemble
IMF is best described by a mean high-mass slope of $\overline{\Gamma} =
1.49\pm0.18$, with an intrinsic scatter of $\sigma^{2}_{\Gamma} =
0.02^{+0.16}_{0.00}$, consistent with a universal IMF. We find no dependence of
the IMF on environmental impacts such as the local star formation rate or
galactocentric radius within M33, which serves as a proxy for metallicity. This
$\overline{\Gamma}$ measurement is consistent with similar measurements in M31,
despite M33 having a much higher star formation rate intensity. While this
measurement is formally consistent with the canonical Kroupa ($\Gamma = 1.30$)
IMF, as well as the Salpeter ($\Gamma = 1.35)$) value, it is the second Local
Group cluster sample to show evidence for a somewhat steeper high-mass IMF
slope. We explore the impacts a steeper IMF slope has on a number of
astronomical sub-fields.
