We derive the time-dependent probability distribution for the number of infected individuals at a given time in a stochastic Susceptible-Infected-Susceptible (SIS) epidemic model. The mean, variance, skewness, and kurtosis of the distribution are obtained as a function of time. We study the effect of noise intensity on the distribution and later derive and analyze the effect of changes in the transmission and recovery rates of the disease. Our analysis reveals that the time-dependent probability density function exists if the basic reproduction number is greater than one. It converges to the Dirac delta function on the long run (entirely concentrated on zero) as the basic reproduction number tends to one from above. The result is applied using published COVID-19 parameters and also applied to analyze the probability distribution of the aggregate number of COVID-19 cases in the United States for the period: January 22, 2020-March 23, 2021. Findings show that the distribution shifts concentration to the right until it concentrates entirely on the carrying infection capacity as the infection growth rate increases or the recovery rate reduces. The disease eradication and disease persistence thresholds are calculated.