We developed an analytic theory of inhomogeneous superconducting pairing in strongly disordered materials, which are moderately close to superconducting-insulator transition. Single-electron eigenstates are assumed to be Anderson-localized, with a large localization volume. Superconductivity develops due to coherent delocalization of originally localized pre-formed Cooper pairs. The key assumption of the theory is that each such pair is coupled to a large number Z of similar neighboring pairs. We derived integral equations for the probability distribution of local superconducting order parameter \Delta and analyzed their solutions in the limit of small dimensionless Cooper coupling constant. The shape of the order parameter distribution is found to depend crucially upon the effective number of "nearest neighbors" defined as Z \nu_0 \Delta. The solution we provide is valid both at large and small effective numbers of neighbors; the latter case is nontrivial as the distribution function is heavily non-Gaussian. Analytic results are supplemented by numerical data, and good agreement between them is observed.