Strength functions have been studied numerically for dense interacting boson systems (with the example of 10 bosons in 5 single-particle states) using a Hamiltonian H, which is a sum of mean-field one-body H1 and a random two-body interaction H2 with strength λ. The strength functions Fξ(E) (defined with respect to H1 basis) changes gradually form Breit–Wigner (BW) to Gaussian form in the chaotic domain as λ increases. A function interpolating these two forms is shown to describe the intermediate region and to determine λF that marks the transition to Gaussian form. Also given is a complete analytical description of the variance of the strength function. The interpolating form for the strength function is used to study the variation of occupancy of single-particle states with excitation energy E.