In the field of Maritime Autonomous Surface Ships (MASS), the accurate modeling of ship maneuvering motion for harbor maneuvers is a crucial technology. Non-parametric system identification (SI) methods, which do not require prior knowledge of the target ship, have the potential to produce accurate maneuvering models using observed data. However, the modeling accuracy significantly depends on the distribution of the available data. To address these issues, we propose a probabilistic prediction method of maneuvering motion that incorporates ensemble learning into a non-parametric SI using feedforward neural networks. This approach captures the epistemic uncertainty caused by insufficient or unevenly distributed data. In this paper, we show the prediction accuracy and uncertainty prediction results for various unknown scenarios, including port navigation, zigzag, turning, and random control maneuvers, assuming that only port navigation data is available. Furthermore, this paper demonstrates the utility of the proposed method as a maneuvering simulator for assessing heading-keeping PD control. As a result, it was confirmed that the proposed method can achieve high accuracy if training data with similar state distributions is provided, and that it can also predict high uncertainty for states that deviate from the training data distribution. In the performance evaluation of PD control, it was confirmed that considering worst-case scenarios reduces the possibility of overestimating performance compared to the true system. Finally, we show the results of applying the proposed method to full-scale ship data, demonstrating its applicability to full-scale ships.