Dark matter structures are expected to exist over a large range of scales, and their properties and distribution can strongly correlate with the underlying particle physics. In this talk, I will describe two separate methods to statistically infer the properties of dark matter substructure using (astrometric)-weak and strong lensing observations, respectively. In the first part of the talk, I will describe how the motion of subhalos in the Milky Way induces a correlated pattern of motions in background celestial objects---known as astrometric weak lensing---and how global signatures of these correlations can be measured using the vector spherical harmonic decomposition formalism. These measurement can be used to statistically infer the nature of substructure, and I will show how this can be practically achieved with future astrometric surveys and/or radio telescopes such as WFIRST and the Square Kilometer Array. Next, I will describe a novel method to disentangle the collective imprint of dark matter substructure on extended arcs in galaxy-galaxy strong lensing systems using likelihood-free (or simulation-based) inference techniques. This method uses neural networks to directly estimate the likelihood ratios associated with population-level parameters characterizing substructure within lensing systems. I will show how this method can provide an efficient and principled way to mine the large sample of strong lenses that will be imaged by future surveys like LSST and Euclid to look for signatures of dark matter substructure. I will emphasize how the statistical inference of substructure using these techniques can be used to stress-test the Cold Dark Matter paradigm and probe alternative scenarios such as scalar field dark matter and enhanced primordial fluctuations.