We examine the nonlinear multiplicative stochastic behavior of chaotic inflation with an emphasis on possible non-Gaussian statistics in initial conditions for cosmological large-scale structure formation. Large-scale mode-mode couplings are analyzed using stochastic dynamics. We reach the following conclusions: 1. Coarse-grained (long-wavelength) scalar fields become nonlinear stochastic variables whose evolution shows behavior unexpected in the classical analysis. The interplay between the classical roll-down (drift) and quantum mechanical fluctuations (diffusion) makes the evolution of the scalar fields extremely nonclassical. Only during the very late stages of chaotic inflation do the scalar fields acquire their classical interpretation (i.e., as deterministic variables). 2. The statistics of initial conditions for cosmological density fluctuations depend on the details of scalar field dynamics during inflation. A non-Gaussian distribution of density fluctuations is a generic feature in chaotic inflation models. However, the astrophysical importance of this effect is strongly model dependent. Generally, deviations from Gaussian statistics are dependent on the strength of the nonlinear self-interaction of the scalar fields. Adiabatic fluctuations were expected to show significant non-Gaussian effects in single shot inflation only on superhorizon scales due to the strong constraint imposed by the cosmic microwave background radiation anisotropy limit. This conclusion is valid for a wide range of models. 3. In a simple chaotic double inflation model, non-Gaussian effects can be astrophysically significant, especially on very large scales. In this model, non-Gaussian statistics can be weakly scale dependent. Non-Gaussian effects could be observable on large scales at a level of approximately a few standard deviations. We also find that some chaotic potentials, designed to give non-scale-invariant density fluctuations, result in strong non-Gaussian phase correlations. 4. In multiple field models, non-Gaussian effects can be important. This depends on the nature of the secondary fields (massive or massless). Multiple field dynamics can provide non-Gaussian statistics as well as nonflat (scale-dependent) spectra. We point out that initial conditions for the large-scale structure in inflationary scenarios in the simplest models are Gaussian to a very high precision. However, as we are forced to consider more finely tuned models of inflation, Gaussian initial conditions should not be taken for granted, even on astrophysically interesting (i.e., observable) scales.