A new approach is proposed for numerical calculations of optical line shapes for disordered molecular aggregates using the excitonic Hamiltonian. The method proposed previously for Gaussian disorder is generalized. A fluctuating origin of the energy scale is introduced as a new independent variable. The integration over this variable eliminates the energy delta function from the integrand, which drastically simplifies the Monte Carlo procedure. The new approach is suitable for calculations using both the diagonalization of the Hamiltonian and the resolvent of the Hamiltonian. As an example, optical line shapes have been calculated for aggregates with bi-Gaussian disorder.