Statical nuclear spectroscopy (also called spectral distribution method), introduced by J.B.French in late 60's and developed in detail in the later years by his group and many other groups, is based on the Gaussian forms for the state (eigenvalue) and transition strength densities in shell model spaces with their extension to partial densities defined over shell model subspaces. The Gaussian forms have their basis in embedded random matrix ensembles with nuclear Hamiltonians consisting of a mean-field one-body part and a residual two-body part. However, following the recent random matrix results for the so called Sachdev-Ye-Kitaev model due to Verbaaarschot et al, embedded random matrix ensembles with k-body interactions are re-examined and it is shown that the density of states, transition strength densities and strength functions (partial densities) in fact follow more closely the q-normal distribution (the parameter q is related to the fourth moment of these distributions with q = 1 giving Gaussian and q = 0 giving semi-circle form). The q-normal has the important property that it is bounded for 0 ≤ q < 1. The q-normal (also its bivariate and general multi-variate extensions) and the associated q-Hermite polynomials are studied for their properties by Bryc, Szabowski and others [P.J. Szabowski, Electronic Journal of Probability 15, 1296Probability 15, (2010]. Following these, in the present article developed is statistical nuclear spectroscopy based on q-normal (univariate and bivariate) distributions and the associated q-Hermite polynomials. In particular, formulation is presented for nuclear level densities, shell model orbit occupancies, transition strengths (for electromagnetic and β and double β-decay type operators) and strength sums.