A harmonic responsibility estimation method based on the Gaussian mixture model (GMM) is proposed for partially observable systems. This method estimates the harmonic responsibility of each harmonic load based on the probabilistic distribution characteristics of measured harmonic voltages, circumventing the difficulties in quantifying harmonic responsibility due to the introduction of unmeasurable line parameters. Specifically, the process begins by training a GMM using the measured harmonic voltage samples. Then, the number and range of Gaussian components in the mixture model are determined based on the Bayesian information criterion and the Kullback-Leibler divergence ratio. Additionally, anomaly detection of harmonic voltage samples is achieved through the Z-test principle. Finally, the effectiveness of the proposed method is verified using the IEEE 14-node test system.