We study the multifractal analysis ͑MFA͒ of electronic wave functions at the localization-delocalization transition in the three-dimensional Anderson model for very large system sizes up to 240 3 . The singularity spectrum f͑␣͒ is numerically obtained using the ensemble average of the scaling law for the generalized inverse participation ratios P q , employing box-size and system-size scaling. The validity of a recently reported symmetry law ͓Mirlin et al., Phys. Rev. Lett. 97, 046803 ͑2006͔͒ for the multifractal spectrum is carefully analyzed at the metal-insulator transition. The results are compared to those obtained using different approaches, in particular the typical average of the scaling law. System-size scaling with ensemble average appears as the most adequate method to carry out the numerical MFA.