This article proposes a three-dimensional (3D) stochastic progressive damage simulation model for the fiber-reinforced polymer matrix composite laminates. Constituents’ strength and fracture energy of each lamina were accounted as spatially varying random fields by 3D Karhunen–Loeve expansion method. Damage mechanisms modeled by continuum damage mechanics include internal damages of matrix and fiber. An anisotropic damage model was used for damage initiation and evolution of each lamina. By the proposed model, we investigated influences of variabilities in the lamina strength properties, and fracture energies along with delamination damage to the damage progression, strength, and post-failure behavior.