Low frequency or long period (LP) earthquakes are a common phenomenon at active volcanoes, and are ubiquitous at persistently active andesitic and dacitic subduction zone volcanoes. At these systems, LP earthquakes provide critical information regarding the state of volcanic unrest, and their occurrence rates are key data on which eruption forecasts are based. Point process modeling of volcanic earthquake occurrence allows potential insights into the underlying physical processes driving unrest, and quantitative, probabilistic, eruption forecasts, for example, through application of the Failure Forecast Method (FFM). However, unlike high-frequency volcano-tectonic (VT) earthquakes, which are typically random or clustered in time, LPs are more commonly quasi-periodic or \textit{anti-clustered}. Consequently, the existing Poisson point process methods used to model occurrence rates of VT earthquakes are unlikely to be optimal for LP data. Here we evaluate the performance of candidate inhomogeneous point process formulations of the FFM for quasi-periodic LP data, based on four different inter-event time distributions: exponential (for Poisson), gamma, inverse Gaussian, and Weibull. Using example LP data recorded before a large explosion at Tungurahua volcano, Ecuador, we examine how well these models explain the observed data, and the quality of retrospective forecasts of eruption time. We use a Markov chain Monte Carlo approach to estimate parameter posterior distributions within a Bayesian framework. Goodness-of-fit is assessed using Quantile-Quantile and Kolmogorov-Smirnov methods, and results are benchmarked against those obtained from idealized synthetic datasets. Inverse Gaussian and gamma models were both found to fit the data well, with the inverse Gaussian model slightly outperforming the gamma model. However, retrospective forecasting analysis shows that the gamma model performs best, with the initial preference for the inverse Gaussian model controlled by catalog incompleteness late in the sequence. The gamma model fits the data significantly better than the Poisson model, and simulations show it produces better forecasts for highly periodic data. Simulations also show that forecast precision increases with the degree of periodicity of the earthquake process using the gamma model, and so should be better for LP earthquakes than VTs. These results provide a new framework for point process modeling of volcanic earthquake time series, and verification of alternative models.