We study the significance of non-Gaussianity in the likelihood of weak lensing shear two-point correlation functions, detecting significantly non-zero skewness and kurtosis in one-dimensional marginal distributions of shear two-point correlation functions in simulated weak lensing data though the full multivariate distributions are relatively more Gaussian. We examine the implications in the context of future surveys, in particular LSST, with derivations of how the non-Gaussianity scales with survey area. We show that there is no significant bias in one-dimensional posteriors of $Ω_m$ and $σ_8$ due to the non-Gaussian likelihood distributions of shear correlations functions using the mock data (100 deg$^2$). We also present a systematic approach to constructing an approximate multivariate likelihood function by decorrelating the data points using principal component analysis (PCA). When using a subset of the PCA components that account for the majority of the cosmological signal as a data vector, the one-dimensional marginal likelihood distributions of those components exhibit less skewness and kurtosis than the original shear correlation functions. We further demonstrate that the difference in cosmological parameter constraints between the multivariate Gaussian likelihood model and more complex non-Gaussian likelihood models would be even smaller for an LSST-like survey due to the area effect. In addition, the PCA approach automatically serves as a data compression method, enabling the retention of the majority of the cosmological information while reducing the dimensionality of the data vector by a factor of ∼5.