A paradigm for constructing and analyzing non-Poisson stimulus-response models of neural spike train activity is presented. Inhomogeneous gamma (IG) and inverse Gaussian (IIG) probability models are constructed by generalizing the derivation of the inhomogeneous Poisson (IP) model from the exponential probability density. The resultant spike train models have Markov dependence. Quantile-quantile (Q-Q) plots and Kolmogorov-Smirnov (K-S) plots are developed based on the rate-rescaling theorem to assess model goodness-of-fit. The analysis also expresses the spike rate function of the neuron directly in terms of its interspike interval (ISI) distribution. The methods are illustrated with an analysis of 34 spike trains from rat CA1 hippocampal pyramidal neurons recorded while the animal executed a behavioral task. The stimulus in these experiments is the animal's position in its environment and the response is the neural spiking activity. For all 34 pyramidal cells, the IG and IIG models gave better fits to the spike trains than the IP. The IG model more accurately described the frequency of longer ISIs, whereas the IIG model gave the best description of the burst frequency, i.e. ISIs < or = 20 ms. The findings suggest that bursts are a significant component of place cell spiking activity even when position and the background variable, theta phase, are taken into account. Unlike the Poisson model, the spatial and temporal rate maps of the IG and IIG models depend directly on the spiking history of the neurons. These rate maps are more physiologically plausible since the interaction between space and time determines local spiking propensity. While this statistical paradigm is being developed to study information encoding by rat hippocampal neurons, the framework should be applicable to stimulus-response experiments performed in other neural systems.