Modern analysis on parton distribution functions (PDFs) requires calculations of the log-likelihood functions from thousands of experimental data points, and scans of multidimensional parameter space with tens of degrees of freedom. In conventional analysis the Hessian method has been widely used to approximate profile of the log-likelihood functions in neighborhood of the best-fit, and to estimate the PDF uncertainties. We propose to use Neural Networks (NNs) and machine learning techniques to model profile of the log-likelihood functions or cross sections for multi-dimensional parameter space, which works beyond the quadratic approximations and meanwhile ensures efficient scans of the full parameter space. We demonstrate efficiency of the new approach in the framework of CT18 global analysis of PDFs by constructing NNs for various target functions, and performing Lagrangian multiplier scans on PDFs and cross sections at hadron colliders. We further study impact of the NOMAD dimuon data on constraining PDFs with the new approach, and find enhanced strange-quark distributions and reduced PDF uncertainties. Moreover, we show how the approach can be used to constrain new physics beyond the Standard Model (BSM) by a joint fit of both PDFs and Wilson coefficients of operators in the SM effective field theory.