Phylodynamics is a set of population genetics tools that aim at reconstructing demographic history of a population based on molecular sequences of individuals sampled from the population of interest. One important task in phylodynamics is to estimate changes in (effective) population size. When applied to infectious disease sequences such estimation of population size trajectories provides information about changes in the number of infections. Current phylodynamic methods either use non-parametric approaches (e.g., curve-fitting based on Gaussian Processes) or model the number of infected individuals deterministically using differential equations. The former methods yield hard-to-interpret results. The latter methods provide estimates of important epidemiological parameters, such as infection and recovery rates, but ignore variation in the dynamics of infectious disease spread. We propose a Bayesian model that combines phylodynamic inference and stochastic epidemic models. In our framework, the population trajectory is modeled by a stochastic compartmental epidemic model, e.g., canonical Susceptible-Infected-Recovered model. To achieve computational tractability, we apply the Linear Noise Approximation approach to approximate stochastic epidemic models with multivariate Gaussian processes (GP). We use modern Markov chain Monte Carlo tools for sampling latent GP models to approximate the posterior distribution of the disease transmission parameters and latent trajectories of compartment sizes. In a simulation study, we show that our method can successfully recover parameters of stochastic epidemic models. We apply our estimation technique to Ebola genealogy estimated using viral genome data from the 2014 Sierra Leone epidemic. We estimate the population trajectory of susceptible and infected individual and changes of the Ebola infection rate through time.