Reconstructing the properties of the astrophysical population of binary compact objects in the universe is a key science goal of gravitational wave detectors. This goal is hindered by the finite strain, frequency sensitivity and observing time of current and future detectors. This implies that we can in general observe only a selected subset of the underlying population, with limited event statistics, and also nontrivial observational uncertainties in the parameters of each event. In this work, we will focus on observations of massive black hole binaries with the Laser Interferometer Space Antenna (LISA). If these black holes grow from population III star remnants ("light seeds"), a significant fraction of the binary population at low masses and high redshift will be beyond LISA's observational reach; thus, selection effects have to be accounted for when reconstructing the underlying population. Here we propose an iterative, kernel density estimation-based nonparametric method, in order to tackle these statistical challenges in reconstructing the astrophysical population distribution from a finite number of observed signals over total mass and redshift. We test the method against a set of simulated LISA observations in a light seed formation scenario. We find that our approach is successful at reconstructing the underlying astrophysical distribution in mass and redshift, except in parameter regions where zero or order(1) signals are observed.
