We estimate the scattering matrix of an arbitrarily complex linear, passive, time-invariant system with $N$ monomodal lumped ports by inputting and outputting waves only via a fixed set of $N_\mathrm{A}<N$ ports while terminating the remaining $N_\mathrm{S}=N-N_\mathrm{A}$ "not-directly-accessible" (NDA) ports with tunable individual loads. First, we present a closed-form approach requiring at least three arbitrary, distinct, and known loads at each NDA port; sign ambiguities on off-diagonal scattering coefficients associated with NDA ports are inevitable. Being matrix-valued, our approach is ideally suited to mitigate noise sensitivity using more accessible ports. It also yields $1+2N_\mathrm{S}+N_S(N_S-1)/2$ as upper bound on the number of required measurements $N_\mathrm{cal}$ for $N_\mathrm{A}>1$ in the low-noise regime. Second, we present a gradient-descent approach using random load configurations, enabling flexible adjustments of $N_\mathrm{cal}$ to further mitigate noise. Third, we present an intensity-only gradient-descent approach, dispensing with phase-sensitive detectors at the expense of an additional blockwise phase ambiguity. Then, we discuss in what applications the inevitable remaining ambiguities are problematic and how to lift them. Finally, we experimentally validate all three approaches with an eight-port reverberation chamber and $N_\mathrm{A}=N_\mathrm{S}=4$, systematically assessing the sensitivity to noise and $N_\mathrm{cal}$. We coin our technique "virtual vector network analyzer (VNA)" because it implies that suitably tunable and characterized individual loads can essentially be interpreted as additional "virtual" VNA ports. Our method can characterize static large antenna systems with many and/or embedded ports, but also reconfigurable wave systems; it may further enable wireless sensing in indoor surveillance, non-destructive testing, and bioelectronics.