We build a statistical framework to infer the global properties of the satellite system of the Andromeda galaxy (M31) from the properties of individual dwarf galaxies located in the Pan-Andromeda Archaelogical Survey (PAndAS) and the previously determined completeness of the survey. Using forward modeling, we infer the slope of the luminosity function of the satellite system, the slope of its spatial density distribution, and the size-luminosity relation followed by the dwarf galaxies. We find that the slope of the luminosity function is β = -1.5 ± 0.1. Combined with the spatial density profile, it implies that, when accounting for survey incompleteness, M31 hosts ${92}_{-26}^{+19}$ dwarf galaxies with M V < -5.5 and a sky-projected distance from M31 between 30 and 300 kpc. We conclude that many faint or distant dwarf galaxies remain to be discovered around Andromeda, especially outside the PAndAS footprint. Finally, we use our model to test if the higher number of satellites situated in the hemisphere facing the Milky Way could be explained simply by the detection limits of dwarf galaxy searches. We rule this out at >99.9% confidence and conclude that this anisotropy is an intrinsic feature of the M31 satellite system. The statistical framework we present here is a powerful tool to robustly constrain the properties of a satellite system and compare those across hosts, especially considering the upcoming start of the Euclid or Rubin large photometric surveys that are expected to uncover a large number of dwarf galaxies in the Local Volume.