The fitting of the observed redshifts and magnitudes of type Ia supernovae to what we would see in homogeneous cosmological models has led to constraints on cosmological parameters. However, in doing such fits it is assumed that the sampled supernovae are moving with the Hubble flow, i.e. that their peculiar velocities are zero. In reality, peculiar velocities will modify supernova data in a way that can impact best-fit cosmological parameters. We theoretically quantify this effect in the nonlinear regime with a Monte-Carlo analysis, using data from semi-analytic galaxy catalogs that are built from the Millennium N-body simulation. We find scaling relations for the errors in best-fit parameters resulting solely from peculiar velocities, as a function of the total number of sources in a supernova survey N and its maximum redshift z_max. For low redshift surveys, we find that these errors can be of the same order of magnitude as the errors due to an intrinsic magnitude scatter of 0.1 mag. For a survey with N=2000 and z_max=1.7, we estimate that the expected peculiar velocity-induced errors in the best-fit cosmological constant density and equation of state can be sigma_Lambda~0.009 and sigma_w~0.01, respectively, which are subdominant to the errors due to the intrinsic scatter. We further find that throwing away supernova data below a redshift z~0.01-0.02 can reduce the combined error, due to peculiar velocities and the intrinsic scatter, but by only about 10%.