The solutions to Burgers equation, in the limit of vanishing viscosity, are investigated when the initial velocity is a Brownian motion (or fractional Brownian motion) function, i.e. a Gaussian process with scaling exponent 0 < h < 1 (type A) or the derivative thereof, with scaling exponent -1 < h < 0 (type B). Large-size numerical experiments are performed, helped by the fact that the solution is essentially obtained by performing a Legendre transform. The main result is obtained for type A and concerns the Lagrangian function x(a) which gives the location at time t = 1 of the fluid particle which started at the location a. It is found to be a complete Devil's staircase. The cumulative probability of Lagrangian shock intervals DELTA-a (also the distribution of shock amplitudes) follows a (DELTA-a)-h law for small DELTA-a. The remaining (regular) Lagrangian locations form a Cantor set of dimension h. In Eulerian coordinates, the shock locations are everywhere dense. The scaling properties of various statistical quantities are also found. Heuristic interpretations are provided for some of these results. Rigorous results for the case of Brownian motion are established in a companion paper by Ya. Sinai. For type B initial velocities (e.g. white noise), there are very few small shocks and shock locations appear to be isolated. Finally, it is shown that there are universality classes of random but smooth (non-scaling) initial velocities such that the long-time large-scale behavior is, after rescaling, the same as for type A or B.