This work emphasizes on predicting probability density function of damages or "number density of damage" in graphite/epoxy polymer matrix composite materials (PMC) under hygrothermal aging condition. A coupled Forward-Backward Stochastic Differential Equation (FBSDE) is proposed as a mathematical model to predict number density of damages. The FBSDE consists of damage nucleation and annihilation rate in terms of Brownian motion. The uncertainty in damage nucleation and annihilation rate is noticed by proposing these two terms as "Brownian motion with drift". In order to verify the proposed model, a quantitative analysis was carried out on a limited number of graphite/epoxy specimens manufactured by VARTM process. The specimens were kept in a hygrothermal condition with room temperature cycling. A rigorous quantitative analysis of damages was done by optical microscopic inspection at different stages of aging period. The damages were classified based on the criteria of their size. Finally, the experimentally collected data for number density of damages were verified by using the proposed FBSDE. A detail parametric study was carried out using FBSDE and best possible predicted data was validated with the experimental observation. A reasonable estimation was observed from the model output.