The quantum spectral statistical properties, i.e. the nearest neighbor level spacing distribution (NNSD), the mode fluctuation distribution (MFD), the spectral rigidity, etc., of a spin-1/2 particle in three-dimensional coupled quartic oscillator potentials are numerically studied. Changing coupling parameters, the system is continuously transformed from an integrable one to chaotic ones. Various kinds of the Gaussian ensembles (GE’s) and the intermediate systems between the integrable class and each GE class can be achieved by selecting non-zero coupling parameters. Especially, this system gives a rare example of the Gaussian symplectic ensemble class. In order to have reliable statistics of quantum levels, it is necessary to evaluate thousands of energy levels from the ground state without missing. We compute the quantum energy levels by numerical diagonalization of the truncated matrix of the Hamiltonian in the basis of harmonic oscillators. Several types of the interpolation formulae of the Poissonian and the Wigner-like distribution of each GE for the NNSD are examined. It is found that the MFD always becomes Gaussian, if the system is chaotic, irrespective of the kinds of the GE’s. Study of statistical properties of quantum levels has started from 1950’s. The random matrix theory (RMT) was first applied to the spectra of neutron resonances in the nuclear scattering experiment. 1) The RMT can be also applicable to predict the distribution of quantum levels in chaotic regime. Nearest neighbor level spacing distribution (NNSD), spectral rigidity, mode fluctuation distribution, distribution of the values of the wave functions, level curvatures with respect to internal parameters, etc., have been shown to be useful for studying the quantum chaos. Among them the NNSD is best known as a measure of the chaoticity and the integrability. It becomes the Wigner distribution if the systems are chaotic, whereas it becomes the Poisson distribution for integrable systems. For intermediate systems between chaotic and integrable regions, we have not reached a satisfactory surmise yet. The several formulae which are based on a physical inspection have been proposed and tested. However, at least in the low energy limit the most successful one is still the Brody distribution. Most of quantum levels which are numerically and experimentally obtained stay in a relatively lower energy region. The Brody distribution is essentially an interpolation formula between the Poissonian and the Wigner surmise and has no particular physical motivation. 2) The Berry-Robnik distribution is a very precious exception with the proper physical inspection. 3) It was introduced under the assumption that the phase space can be separated into two regions, that is, the integrable region (the KAM tori area) and