This paper introduces a novel framework to construct the region of attraction (ROA) of a power system centered around a stable equilibrium by using stable state trajectories of system dynamics. Most existing works on estimating ROA rely on analytical Lyapunov functions, which are subject to two limitations: the analytic Lyapunov functions may not be always readily available, and the resulting ROA may be overly conservative. This work overcomes these two limitations by leveraging the converse Lyapunov theorem in control theory to eliminate the need of an analytic Lyapunov function and learning the unknown Lyapunov function with the Gaussian Process (GP) approach. In addition, a Gaussian Process Upper Confidence Bound (GP-UCB) based sampling algorithm is designed to reconcile the trade-off between the exploitation for enlarging the ROA and the exploration for reducing the uncertainty of sampling region. Within the constructed ROA, it is guaranteed in probability that the system state will converge to the stable equilibrium with a confidence level. Numerical simulations are also conducted to validate the assessment approach for the ROA of the single machine infinite bus system and the New England $39$-bus system. Numerical results demonstrate that our approach can significantly enlarge the estimated ROA compared to that of the analytic Lyapunov counterpart.