Nonlinear hyperbolic balance laws with uncertain (random) initial data are ubiquitous in a plethora of transport phenomena that often exhibit shocks. We develop the method of distributions for such problems by adding a model error term to a deterministic equation for the cumulative distribution function (CDF) of the system states. We use two alternative strategies, Newtonian relaxation and neural networks, to infer this term from observations of the system dynamics. The former strategy is amenable to theoretical analysis of its convergence with respect to data sparsity, while the latter offers more flexibility. The CDF equation is exact for linear conservation laws and nonlinear conservation laws with a smooth solution, such that the CDF equation can be used to formulate predictions at times when observations cease to be available. Whenever shocks develop as a result of the nonlinearity, observations are used to detect the discrepancy that emerges as model error. Spatial data density is crucial for good interpolation accuracy, whereas long temporal sequences of observations improve future projections.