  Simulations of Type Ia Supernovae (SNIa) surveys are a critical tool for
correcting biases in the analysis of SNIa to infer cosmological parameters.
Large scale Monte Carlo simulations include a thorough treatment of observation
history, measurement noise, intrinsic scatter models and selection effects. In
this paper, we improve simulations with a robust technique to evaluate the
underlying populations of SNIa color and stretch that correlate with
luminosity. In typical analyses, the standardized SNIa brightness is determined
from linear `Tripp' relations between the light curve color and luminosity and
between stretch and luminosity. However, this solution produces Hubble residual
biases because intrinsic scatter and measurement noise result in measured color
and stretch values that do not follow the Tripp relation. We find a $10\sigma$
bias (up to 0.3 mag) in Hubble residuals versus color and $5\sigma$ bias (up to
0.2 mag) in Hubble residuals versus stretch in a joint sample of 920
spectroscopically confirmed SNIa from PS1, SNLS, SDSS and several low-z
surveys. After we determine the underlying color and stretch distributions, we
use simulations to predict and correct the biases in the data. We show that
removing these biases has a small impact on the low-z sample, but reduces the
intrinsic scatter $\sigma_{\textrm{int}}$ from $0.101$ to $0.083$ in the
combined PS1, SNLS and SDSS sample. Past estimates of the underlying
populations were too broad, leading to a small bias in the equation-of-state of
dark energy $w$ of $\Delta w=0.005$.
