In this paper, we present a unified mathematical framework for neural symbolic (NeSy) systems, called Neural Symbolic Energy-Based Models (NeSy-EBMs). NeSy-EBMs encompass both discriminative and generative NeSy modeling, and provide general formulas for gradient estimation of key learning losses. We present four training methods that leverage methods from various fields, including bi-level and probabilistic policy optimization. We implement the NeSy-EBM framework with Neural Probabilistic Soft Logic (NeuPSL), an open-source NeSy-EBM library designed for scalability and expressiveness, and demonstrate the practical benefits of NeSy-EBMs for various tasks, including image classification, graph node labeling, autonomous vehicle situational awareness, and question answering.