this paper, we first introduce a symmetric phase type (PH) distribution, and then propose a symmetric exponential-polynomial type (EP) distribution based on the symmetric PH distribution. The class of symmetric EP distributions is illustrated to be so large that an arbitrary symmetric probability density function in L2 (-∞, +∞), which is the space of square integrable functions on the real line, can be approximated in (-∞, +∞) by a sequence of symmetric EP probability density functions. We prove this result by means of a constructive approach based on two orthogonal decompositions: Laguerre spectrum decomposition and Hermite spectrum decomposition. We also provide a moment-based approach for simply determining a symmetric EP probability density function to approximatively express a symmetric random variable under which some moments of the symmetric random variable are given. We further propose multidimensional symmetric PH and EP distributions and provide their useful structures and properties. Finally, we apply the symmetric PH and EP distributions to study the probability Hough transform.