In a network, the shortest paths between nodes are of great importance as they allow the fastest and strongest interaction between nodes. However measuring the shortest paths between all nodes in a large network is computationally expensive. In this paper we propose a method to estimate the shortest path length (SPL) distribution of a network by random walk sampling. To deal with the unequal inclusion probabilities of dyads (pairs of nodes) in the sample, we generalize the usage of Hansen-Hurwitz estimator and Horvitz-Thompson estimator (and their ratio forms) and apply them to the sampled dyads. Based on theory of Markov chains we prove that the selection probability of a dyad is proportional to the product of the degrees of the two nodes. To approximate the actual SPL for a dyad, we use the observed SPL in the induced subgraph for networks with large degree variability, i.e., the standard deviation is at least two times of the mean, and for networks with small degree variability, estimate the SPL using landmarks for networks with small degree variability. By simulation studies and applications to real networks, we find that 1) for large networks, high estimation accuracy can be achieved by using a single random or multiple random walks with total number of steps equal to at least 20% of the nodes in the network; 2) the estimation performance increases as the network size increases but tends to stabilize when the network is large enough; 3) a single random walk performs as well as multiple random walks; 4) the Horvitz-Thompson ratio estimator performs best among the four estimators.