Bosonic quantum fields and quantum devices can be simulated in phase space, to overcomes exponential complexity issues as the Hilbert space dimensionality grows. Algorithms for sampling initial quantum states are analyzed here as a fundamental issue in these simulations. The generalized P, Wigner, and Q-function phasespace methods are treated, with applications ranging from Bose-Einstein condensates to quantum computers and photonic networks. Quantum initial conditions treated include multimode Gaussian states and general number state expansions. These are sampled probabilistically, in complex-P-distribution cases with additional complex weights. For high-order correlations of the type found in many quantum technology experiments, different representations and sampling techniques have a sampling error that changes exponentially as the mode number increases. Purely probabilistic techniques may scale exponentially worse than experimental shot-noise error scaling. Yet, combining random samples with complex weighting techniques gives error scaling that is equal or even exponentially better than experimental shot-noise scaling. We show that complex weighted P-distribution sampling can also be used to initialize simulations with the Wigner distribution, which is difficult to sample owing to its nonpositive character. As a result, complex weighted sampling can require exponentially less samples than either experiment or unweighted probabilistic sampling, leading to efficient simulation techniques for quantum networks and BECs.