{'fraction_non_alphanumeric': 0.040537642096232215, 'fraction_numerical': 0.01508543647672839, 'mean_word_length': 4.421705426356589, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Sparse spectral estimation from point process observations

We consider the problem of estimating the power spectral density of the neural covariates underlying the spiking of a neuronal population. We assume the spiking of the neuronal ensemble to be described by Bernoulli statistics. Furthermore, we consider the conditional intensity function to be the logistic map of a second-order stationary process with sparse frequency content. Using the binary spiking data recorded from the population, we calculate the maximum a posteriori estimate of the power spectral density of the process while enforcing sparsity-promoting priors on the estimate. Using both simulated and clinically recorded data, we show that our method outperforms the existing methods for extracting a frequency domain representation from the spiking data of a neuronal population.


INTRODUCTION
Spectral analysis of electroencephalography (EEG) recordings has long been used for research and diagnosis purposes such as the identification of sleep disorders [1,2] and epileptic seizures [3,4]. With the development of invasive recording procedures, single-and multi-unit recordings have increasingly become popular in neural sciences [5,6]. These data do not suffer from the low spatial resolution of non-invasive measurements such as EEG, as they capture the spiking activity of a localized population neurons; however, their binary nature introduces various signal processing challenges [7].
Recently, the theory of point processes has provided a mathematical framework to model and analyze neuronal data such as single and multi unit recordings in the time domain [8,9,10]. Following the common frequency-domain analysis of neural data such as EEG, existing methods for point process spectral estimation often compute a continuous estimate of the spiking rate and analyze the power spectral density (PSD) of this estimate. The spiking rate estimation is either done by simply smoothing the spiking histogram [11,12,13] or using generalized linear Gaussian state-space models to estimate the conditional intensity function (CIF) of the point process [14,15]. However, these approaches suffer from the following shortcomings: Firstly, they are limited in terms of their spectral resolution as smoothing in the time domain for spiking rate estimation results in distortion in the frequency domain. Secondly, second-order statistics of the underlying time-series are needed for spectral estimation, and even if the spiking rate is estimated correctly, the second-order statistics need not be. Thirdly, current methods do not account for spectral sparsity which is often observed in the frequency domain analysis of biological data such as EEG.
In this work, we address the above shortcomings by casting the point process spectral estimation problem as a discrete-parameter harmonic spectral estimation in which the the second moments of the harmonic process driving the point process are estimated. To this end, we consider a conditional Bernoulli point process framework and model its CIF as the logistic map of a harmonic process. Then, we use an Expectation-Maximization (EM) algorithm to compute maximum a posteriori (MAP) estimate of the PSD of the harmonic process as the spectral representation of point process using the observed spiking data and sparsity-promoting priors. Simulation results and application to real multi-unit data recorded under general anesthesia verifies the superior performance of our method compared to existing techniques.
The rest of the paper is organized as follows: Section 2 introduces our model for the spiking activity of a population of neurons driven by a harmonic process. In Section 3, we derive the sparse MAP estimator of the PSD associated with the harmonic process. Section 4 provides results for simulated data and real multi-unit recordings under general anesthesia, comparing our sparse PSD estimates with existing methods. This is followed by our concluding remarks in Section 5.

PROBLEM FORMULATION
A point process can be fully characterized by its Conditional Intensity Function (CIF). For t ∈ (0, T ], the CIF of a point process is defined as [16] λ(t|H t ) := lim where H t and N (t) denote respectively the spiking history in (0, t) and the number of spikes in (0, t]. If the point process is discretized by small enough bins of size ∆, the data in each bin (n k ) can be approximated by a Bernoulli random variable with the success probability of λ k := λ k∆|H k∆ ∆, for 0 ≤ k ≤ K, where K = T ∆ . Due to the absolute refractory period of neurons, the sampling period of ∆ ∼ 1 ms is typically sufficient to ensure that at most one spike occurs in any bin [8].
The oscillatory patterns in neuronal spiking can be directly attributed to its CIF. Thus, a spectral analysis of the CIF would provide an informative frequency domain representation of the point process. We consider an ensemble of L neurons, whose spiking data is denoted by k denoting the data of neuron in bin k, and assume that they are driven by the same CIF. This assumption is plausible for multi-unit recordings as they are obtained from neighboring neurons [5]. To associate a power spectral density with the mutual CIF, we consider the CIF to be the logistic map of a second-order stationary signal with mean µ. Denoting the samples of the second-order stationary process by x k for 0 ≤ k ≤ K, the complete discretized model can be expressed as: The objective is then to calculate a sparse estimate of the PSD of x k as the spectral representation of the neural spiking data having observed D.

BAYESIAN ESTIMATION OF THE PSD
According to the Spectral Representation Theorem [17], for the zero-mean and second-order stationary signal x k − µ with PSD S(ω), there exists a continuous, orthogonal increment, and complex process Z(ω) such that where the integral is in the Riemann-Stieltjes sense and We consider a discrete approximation to the PSD S(ω) which results in Z(ω) being constant over intervals of length π N for large enough N [17]. Thus, Z(ω) can be replaced by a jump process in [0, π) with jumps of π N (a i + jb i ) at ω i = iπ N for i = 1, 2, . . . , N −1, where a i and b i are random variables and N controls the degree of approximation. Following the Spectral Representation Theorem and the discrete PSD approximation, the PSD at ω i would be S(ω i ) = π 2 N 2 E{a 2 i +b 2 i }.
Furthermore, considering x k to be real, Eq. (3) can be simplified as: we can write Eq. (4) as x = Av. We further assume that the process Z(ω) is Gaussian. Therefore, based on the definition of random vector v, we have v i ∼ N (0, σ 2 i ), for i = 2, . . . , 2N − 1, and due to the orthogonality of the increments of process Z(ω), v i 's are independent. This re- The estimation of µ is not of particular importance as our goal is to infer the oscillatory patterns in the CIF. However, in order to have a consistent prior on all the elements of v, we assume an independent Gaussian prior on µ such that v 1 ∼ N (0, σ 2 1 ). Therefore, an estimate of the parameter vector θ := [σ 2 1 , σ 2 2 , ..., σ 2 2N −1 ] T results in an estimate of the approximated PSD of x k .
To obtain an sparse estimate of θ in the Bayesian estimation setting, we enforce independent exponential priors with parameter γ on the elements of θ, i.e. log f θ (θ) = (2N − 1) log γ −γ where we estimate γ through two-fold cross validation [18]. This log-prior is equivalent to the 1norm penalty of θ which is known to promote sparsity. Thus, the maximum a posterior (MAP) estimate of θ is defined as: As P (D|θ) is an intractable function of θ and the vector v acts as a latent variable in the Bayesian estimation of θ, we employ the Expectation-Maximization (EM) algorithm to calculate the MAP estimate in Eq. (6) [19,20].
The E Step: Denote the parameter vector estimate at the r th iteration of the EM algorithm by θ (r) = σ 2 1 (r) , σ 2 2 (r) , · · · , σ 2 2N −1 (r) T . In the E step, we calculate the expectation of the complete data log-likelihood as: where we have used the conditional independence relation P (D|v, θ) = P (D|v). In (7), the term cnst. represents all terms which are not functions of θ. Calculating the distribution of v|D, θ (r) involves intractable integrals and Monte Carlo methods are not efficient for calculating the term E v|D, θ (r) v 2 i in (7) as they have to be performed in every EM iteration. However, the density of v|D, θ (r) can be approximated by a multivariate Gaussian density N µ [8,10]. Similar to [8,10], we calculate µ (r) v as the mode of the log-density log f v|D, θ (r) (v|D, θ (r) ) in (8) given by: and evaluate the Hessian of this log-density at the calculated [19]. The maximization problem in (8) is concave, and the Hessian is negative definite due to the diagonal terms of the form for k = 1, 2, · · · , K, which is often referred to as the Peristimulus Time Histogram (PSTH), includes all the spiking information we need to calculate µ (r) v due to our mutual CIF assumption. The multivariate Gaussian estimate results in i 's, we maximize Q(θ| θ (r) ) in (7) with respect to θ in the M step. Q(θ| θ (r) ) is a quasi-concave function over the positive orthant with a unique maximizer which is the non-negative solution of a quadratic equation. Hence, the update rule of the EM algorithm becomes:

APPLICATION TO SIMULATED AND REAL DATA
As for our simulated data, we consider the dual-tone signal of the form x(t) = 1.48 cos(2πf 0 t) + 0.685 cos(2πf 1 t) + 0.17n(t) − 5.7 (10) with f 0 = 1 Hz and f 1 = 10 Hz sampled at f s = 300 Hz as the underlying harmonic process in our model, and n(t) representing a zero-mean white Gaussian noise with unit variance. A total of K = 1000 samples are considered for L = 10 neurons. The bias term of −5.7 is chosen to make sure that the resulting spiking rate of 0.056 in PSTH is low enough and consistent with real-world neuronal spiking rates. Figure  1 shows the harmonic process as well as the resulting raster plot of the ensemble according to our model.   Figure 2-(c) illustrates the normalized multitaper estimate [21] of the PSD corresponding to the smoothed PSTH for both the narrow and wide kernels. The multitaper estimate is among the most reliable nonparametric spectral estimation methods. The foregoing method of estimating the PSD of the smoothed PSTH, which we denote by PSTH-PSD, is one of the most commonly-used methods for extracting the spectral representation of the spiking data [11,12,13]. We have considered a frequency spacing of 0.125Hz for the PSD estimates here, which corresponds to N = 1200 in our proposed method, and DC components of the normalized PSD's are eliminated in the figures. The PSD estimate corresponding to the narrow smoothing kernel (Figure 2-(c), green trace) detects the two peaks at 1 Hz and 10 Hz, but has a high variability in higher frequencies. Also, two spurious peaks at 2 Hz and 9 Hz are detected in the PSD. On the other hand, the estimate corresponding to the wide smoothing kernel ( Figure  2-(c), dotted red trace) has a smaller variability but misses the 10 Hz peak. These results show the high sensitivity of the PSTH-PSD approach to the choice of the smoothing kernel and the absence of a criterion for choosing one. The result of our proposed method in Figure 2-(d) with γ opt = 10 −4 (obtained by two-fold cross-validation) fully recovers the two frequency components after 130 EM iterations.
Finally, we apply our algorithm to multi-unit recordings in [22] from a human subject under Propofol-induced general anesthesia. The original sampling rate of 1 kHz was reduced by the factor of 40 to reduce computational complexity, and a time window of 50s was considered (K = 1250). Out of the 41 neurons in the data set, we have considered the ones with at least two spikes in the 50s time frame (L = 27).  Figure 4 shows the PSD estimates using our proposed method as well as the PSTH-PSD estimates. All PSD's have the frequency spacing of 0.02Hz corresponding to N = 625 in our method. Similar to the foregoing simulation results, Figure 4-(c) shows the considerable variability of the narrow kernel PSTH-PSD in higher frequencies while the wide kernel PSTH-PSD suppresses the PSD components above 0.4Hz, implying the high sensitivity of the estimates to the choice of the kernel width. However, our method in Fig. 4-(d), after 100 EM iterations, calculates a sparse PSD with the regularization parameter obtained by two-fold cross-validation. It is worth noting that the results of the state-space model of [10,14,15] were very similar to the wide Gaussian kernel PSTH-PSD [19], and thus have not been shown in Figures 2  and 4 for brevity.

CONCLUSION
Although the theory of point process provides a concrete mathematical framework for analyzing binary data in the time domain, there has not been much work on extending this theory to incorporate a spectral representation for binary data. Existing approaches calculate an estimate of the spiking rate or the CIF driving the point process by smoothing the PSTH in the time domain or fitting a state-space model to the data. These methods have drawbacks such as limiting the resolution in frequency domain due to smoothing in the time domain and not taking into account the underlying sparse spectral structure often observed in biological signals.
In this paper, we have addressed these issues by introducing a model for spectral estimation from spiking data which directly estimates the PSD as the second-order moments of the underlying process from the observed binary data and employs sparsity-promoting priors for estimation. Application to simulated and real data reveals that our proposed method outperforms existing methods for calculating a point process spectral representation. Finally, although we have only focused on the application of our spectral estimation method to neuronal multi-unit recordings, our techniques can be similarly applied to any binary data exhibiting oscillatory behavior.