# Parameter Estimation of a Class of Neural Systems with Limit Cycles
## 4 Parameter Estimation of the Spiking Neurons  
### Least-Squares Estimation Algorithms  
To estimate the parameters \(\theta\), let us consider the cost function  
\[\mathbb{J}(\theta)=\sum_{k=1}^{L_{d}}\sum_{i=1}^{2}(y_{i}(k)-\phi_{i}^{\top}( k)\theta)^{2},\]  
where \(L_{d}\) denotes the data length. Using the least-squares principle to solve the optimization problem, we explore the following recursive least-squares (RLS) parameter estimation formulas [3, 27] to estimate the parameter vector \(\theta(k)\). \[\left\{ \begin{array}{rcl}\hat{\theta}(k)&=&\hat{\theta}(k-1)+L(k)\varepsilon(k), \\ \varepsilon(k)&=&\begin{bmatrix}y_{1}(k)-\phi_{1}^{\top}(k)\hat{\theta}(k-1)\\ y_{2}(k)-\phi_{2}^{\top}(k)\hat{\theta}(k-1)\end{bmatrix},\\ L(k)&=&P(k-1)\phi(k)\left[\lambda+\phi^{\top}(k)P(k-1)\phi(k)\right]^{-1},\\ P(k)&=&(I-L(k)\phi^{\top}(k))P(k-1),\end{array} \right. \tag{6}\]  
where \(\lambda\in(0,1]\) is the forgotten factor, \(P(k)\in\mathbb{R}^{h\times h}\) is the vector of adjustment gains, \(\tilde{n}=6\), \(P(0)=p_{0}I\), and \(p_{0}\) is a large positive number. To initialize the RLS algorithm, the initial value \(\hat{\theta}(0)\) is generally taken to be a zero vector or a small real vector, e.g., \(\hat{\theta}(0)=10^{-6}\mathbf{1}_{h}\) with \(\mathbf{1}_{h}\) being an \(\tilde{n}\)-dimensional column vector whose elements are \(1\). Now, the RLS parameter estimation algorithm for the spiking neuron model is summarized in Algorithm 1. 1. Discretize the FHN model based on the explicit forward Euler method with step size \(T\). 2. Initialization: set the forgotten factor \(\lambda\in(0,1]\) and let \(k=1\), \(\hat{\theta}(0)=\mathbf{1}_{h}/p_{0}\), \(P(0)=p_{0}I\), \(p_{0}=10^{6}\). 3. Collect the measurement state data and determine a data length \(L_{d}\). Then, form the output data \(y(k)\) and the information vector \(\phi(k)\). 4. Compute \(L(k)\) and \(P(k)\) by \[L(k)=\frac{P(k-1)\phi(k)}{\lambda+\phi^{\top}(k)P(k-1)\phi(k)},P(k)=(I-L(k) \phi^{\top}(k))P(k-1),\] respectively. 5. Update the estimate \(\hat{\theta}(k)\) by \(\hat{\theta}(k)=\hat{\theta}(k-1)+L(k)\varepsilon(k)\), where \[\varepsilon(k)= \begin{bmatrix}y_{1}(k)-\phi_{1}^{\top}(k)\hat{\theta}(k-1)\\ y_{2}(k)-\phi_{2}^{\top}(k)\hat{\theta}(k-1)\end{bmatrix} .\]
6. If \(k<L_{d}\), increase \(k\) by \(1\) and go to step (3); otherwise, stop the procedure and output the estimate \(\hat{\theta}(L_{d})\) of the parameter vector \(\theta\). **Algorithm 1** RLS algorithm  
To obtain better estimation accuracy, we modify the algorithm (6) and apply a multiinnovation recursive least-squares (MIRLS) parameter estimation algorithm [7] for (5) by introducing an innovation length \(p\). By defining the information matrix \(\Phi(p,k)\) and stacked output vector \(Y(p,k)\) as  
\[Y(p,k)=[y^{\top}(k),y^{\top}(k-1),\cdots,y^{\top}(k-p+1)]^{\top}\in\mathbb{R} ^{2p}, \tag{7}\]  
\[\Phi(p,k)=[\phi^{\top}(k),\phi^{\top}(k-1),\cdots,\phi^{\top}(k-p+1)]^{\top} \in\mathbb{R}^{6\times 2p}, \tag{8}\]  
the innovation vector \(\mathcal{E}(p,k)\) can be expressed as  
\[\mathcal{E}(k)=Y(p,k)-\Phi^{\top}(p,k)\hat{\theta}(k-1).\]  
From here, we present the following MIRLS iterative formulas  
\[\left\{ \begin{array}{rcl}\hat{\theta}(k)&=&\hat{\theta}(k-1)+L(k)\mathcal{E} (k),\\ \mathcal{E}(k)&=&Y(p,k)-\Phi^{\top}(p,k)\hat{\theta}(k-1),\\ L(k)&=&P(k-1)\Phi(p,k)\Big{[}\lambda+\Phi^{\top}(p,k)P(k-1)\Phi(p,k)\Big{]}^{- 1},\\ P(k)&=&(I-L(k)\Phi^{\top}(p,k))P(k-1).\end{array} \right. \tag{9}\]The MIRLS parameter estimation algorithm for the FHN model is summarized in Algorithm 2. 1. Discretize the FHN model based on the explicit forward Euler method with step size \(T\). 2. Initialization: set the forgotten factor \(\lambda\in(0,1]\) and let \(k=1\), \(\hat{\theta}(0)=\mathbf{1}_{h}/p_{0}\), \(P(0)=p_{0}l\), \(p_{0}=10^{6}\). 3. Collect the measurement state data and determine a data length \(L_{d}\). Then, form the output data \(y(k)\) and the information vector \(\phi(k)\). 4. Given an innovation length \(p\), form \(Y(p,k)\) by (7) and \(\Phi(p,k)\) by (8). 5.