We present a large-scale simulation of the three-dimensional Ising spin glass with Gaussian disorder to low temperatures and large sizes using optimized population annealing Monte Carlo. Our primary focus is investigating the number of pure states regarding a controversial statistic, characterizing the fraction of centrally peaked disorder instances, of the overlap function order parameter. We observe that this statistic is subtly and sensitively influenced by the slight fluctuations of the integrated central weight of the disorder-averaged overlap function, making the asymptotic growth behaviour very difficult to identify. Modified statistics effectively reducing this correlation are studied and essentially monotonic growth trends are obtained. The effect of temperature is also studied, finding a larger growth rate at a higher temperature. Our state-of-the-art simulation and variance reduction data analysis suggest that the many pure state picture is most likely and coherent.