What is the best solution to replace NaN values?

You are right in saying that replacing with a simple mean, mode... is a common but unreliable imputation strategy in many cases. You have in scikit learn some utilities for imputation of missing values (have a look at https://scikit-learn.org/stable/modules/classes.html#module-sklearn.impute) using for instance the knn imputer as an additional strategy.

Take into account you cannot assume your feature of interest follows a normal distribution, so instead you can actually apply a kernel density estimator to model such distribution, see here: http://scikit-learn.org/stable/modules/density.html