This paper is devoted to the numerical analysis of weak Galerkin mixed finite element method (WGMFEM) for solving a heat equation with random initial condition. To set up the finite element spaces, we choose piecewise continuous polynomial functions of degree j+1 with j≥0 for the primary variables and piecewise discontinuous vector-valued polynomial functions of degree j for the flux ones. We further establish the stability analysis of both semidiscrete and fully discrete WGMFE schemes. In addition, we prove the optimal order convergence estimates in L2 norm for scalar solutions and triple-bar norm for vector solutions and statistical variance-type convergence estimates. Ultimately, we provide a few numerical experiments to illustrate the efficiency of the proposed schemes and theoretical analysis.