We proposed a novel Monte Carlo algorithm to sample the output distribution of a neural network. The algorithm combines the Wang-Landau algorithm with Gibbs sampling, in which the Monte Carlo proposals are derived from the gradient information of the neural network. Our Wang-Landau algorithm with “gradient updates” is able to traverse the output space of the neural network more effectively than the original Wang-Landau algorithm with “random updates”. This enables a comprehensive study of the input–output relationship in a broad output space, and provides insights into the robustness of the neural network when generalized to unseen data. We applied the algorithm to study several common neural network models for a binary classification task. Our preliminary investigation indicates that a generative model might be able to learn better features to distinguish the classes and map them to a small, unique output range – an indicator of better model quality.