A numerically efficient method of shooting and bouncing rays for sparse particles (SBR-SP) is presented to analyze a large two-dimensional (2D) array of magnetodielectric circular cylinders which are distributed with any positions, sizes, and media on the xy-plane. Iterative- and closed-forms of the SBR-SP are used to understand the convergence and accuracy characteristics of a large array of cylinders and we also compare them with the results obtained by a mode-matching technique, thus resulting in favorable agreement. Especially, the closed-form of the SBR-SP with GPU acceleration is fast and accurate enough to obtain various echowidth behaviors of a 80×80 square array within 71.8 s. In addition, a statistical analysis for a randomly distributed array with 6,400 elements is conducted in terms of mean and standard deviation and we obtain that initial shooting rays excited by an incident field are strongly dominant in statistical behaviors of scattering patterns. Our algorithm based on the SBR can be extended to three-dimensional (3D) problems where we model a large array of lossy dielectric spheres for the estimation of rain attenuation and signal variation.