# Hamiltonian Monte Carlo with Random Effect for Analyzing Cyclist Crash Severity
## 2 Method
### Hamiltonian Monte Carlo
#### 2.2.2 Model Parameters' Estimations  
In this section, we present the simplified model implementation to estimate parameters, being based on the work in the literature review [21]. The process of parameters' estimation could be divided into 5 main sections and related subsections:  
* General data preparation Preparation of the matrix of input by QR decomposition, orthogonal matrix of Q and upper triangular matrix of R, to improve efficiency of HMC sampling. Creating value of \(\varepsilon\) for all parameters to update the model parameters estimates (e.g., Equation (8)). Auxiliary log density function of momentum would be created by assigning the mean of zeroes and variance-covariance diagonal identity matrix of \(p\times p\), where \(p\) is the number of parameters to be estimated. Non-informative initial values of 0 would be given to all parameters.
* Leapfrog algorithm
* Use of gradients Derivative of full posterior distribution relative to \(\beta\), \(\tau\), and \(\xi\) would be obtained, and saved in a matrix for updating various parameters to be estimated. Based on Equation (7), \(m\) or momentum would be updated, and \(\theta\) would be updated based on Equation (8). It should be noted that the \(\varepsilon\) would be set as small values for all parameters. The leapfrog would be conducted for 10 iterations nested inside the iteration of HMC sampling. The gradient in leapfrog would direct the sampling to the optimum values.
* Updating \(\theta\) and \(m\) In each iteration of leapfrog, the previous value of \(\theta\) and \(m\) would be used and updated.
* Estimating log posterior We are still in the main loop of the first HMC sampling iteration. For this section, values would be estimated for Equation (9), related to its numerator and denominator. The numerator includes the new value of \(\theta\) that is just estimated, while the equation's denominator would use the old value of \(\theta\), which for the first iteration is preassigned values of 0.  
Here, full posterior, as discussed before, is dependent on Equation (18), and other priors' value. For Equation (18), parameters related to random effect, \(u\), need to be prepared. That is based on Equation (12), which itself would be written based on Equation (13) or (14). Now, \(log\ P(\beta|\sigma_{\beta}^{2})\), \(log\ P(\tau)\), \(log\ P(\xi)\), related to Equations (15)-(17), respectively, would be calculated and added to Equation (18). Recall that the process for numerator and denominator would be estimated separately for \(\theta\) and \(\theta_{t-1}\).  
* Estimation of \(\alpha\): Whether to reject or accept the new proposal  
As discussed, one of the advantages of the leapfrog algorithm is to increase the rate of acceptance, compared with random walk. Now, based on Equation (9), in case of \(\alpha>1\), we accept the new proposed \(\theta\) and \(-m\), otherwise we reject it and keep \(\theta_{t-1}\) and \(m_{t-1}\). The values would be saved and used for the next iteration. The process would be iterated for 6000 HMC sampling iterations, each with 10 leapfrog iterations. The final \(\theta\) would bind all \(\theta\)s for further analysis.  
It should be noted that, for the HMC, we flip the distribution by getting the log of probability density. The idea is that we follow the ledge of the distribution, which tends to go to the bottom of the distribution based on the gravity. Therefore, by generating more samples from the bottom of the distribution, in reality, we draw from the pick of the distribution.  
* Results preparation  
The main output of the model is the \(\theta\)s value during iterations, Equation (10). Various quantiles of sampled \(\theta\) would be used after discarding the first N number of burning samples, highlighting certainty or uncertainty associated with parameters' estimates.