This paper primarily establishes an asymptotic variance estimate for smooth linear statistics associated with zero sets of systems of random holomorphic sections in a sequence of positive Hermitian holomorphic line bundles on a compact Kähler manifold $(X, ω)$ in a general non-Gaussian setting. Using this variance estimate and the expected distribution, we derive an equidistribution result for zeros of these random systems, which proves that the smooth positive closed form $ω^{k}$ can be approximated by currents of integration along analytic subsets of $X$ of codimension $k$, $k \in \{1, \ldots, n\}$. The probability measures taken into consideration in this paper are sufficiently general to include a wide range of the measures commonly encountered in the literature, for which we give equidistribution results at the end, such as the standard Gaussian measure, Fubini-Study measure, the area measure of spheres, probability measures whose distributions have bounded densities with logarithmic decaying tails and locally moderate measures among others.