# Hamiltonian Monte Carlo with Random Effect for Analyzing Cyclist Crash Severity
## 2 Method
### Hamiltonian Monte Carlo  
One of the shortcomings of application of random walk in algorithms, such as Metropolis Hastings, is the samplers explore the target distribution extremely slowly [18]. In other words, the Metropolis Hastings (MH) takes a long time with large number of iterations to get to the true posterior density (e.g., often the acceptance rate is less than a quarter of all proposals) [19]. Hamiltonian Monte Carlo (HMC), on the other hand, could be employed by utilizing methods such as leapfrog to generate transition and eliminating the random walk behavior. The HMC algorithm has been identified as a powerful tool, which could efficiently compute model parameters. In the HMC method, we assume that the system has an energy level, \(E\), being dependent on \(\theta\) and \(M\), where \(\theta\) is the model's parameter to be estimated, and \(M\) is the momentum. Energy itself depends on the position of a particle and the particle's kinetic energy. We would solve for the particle with potential energy of \(u\) (\(\theta\)), and kinetic energy of \(k(M)\), where \(u\) (\(\theta\)) is \(-\)log\(f\)(\(\theta\)), and the kinetic energy is 0.5 \(p^{T}M^{-1}p\). Here \(M\), the auxiliary momentum variables, is added, which is known. Therefore, we would have:  
\[H(\theta,m)=u+k=-logf(\theta)+0.5m^{T}M^{-1}m \tag{6}\]  
where \(m\)\(\sim N(0,M)\). As we cannot exactly solve for part of a particle, we end up discretizing the path, so the algorithm would be small enough to make a precise enough approximation. Leapfrog would be used to approximate solutions for Hamiltonian equations by using a step size of \(\varepsilon\). The output of leapfrog is the proposed values of \(\theta^{*}\)and \(m^{*}\). The gradient of the above, with respect to \(m\) and \(\theta\), would be the Hamiltonian equations, which would be solved by the leapfrog algorithm as:  
\[m(t+\varepsilon/2)=m(t)+\varepsilon/2\frac{dm}{dt}=m(t)+\varepsilon/2\
abla_ {\theta}logf(\theta) \tag{7}\]  
\[\theta(t+\varepsilon)=\theta(t)+\varepsilon\frac{d\theta}{dt}=\theta(t)+ \varepsilon M^{-1}m(t+\varepsilon) \tag{8}\]  
where \(\
abla_{\theta}logf(\theta)\) is the gradient of log posterior density. Here, few steps would be used to move to a next proposal. In summary, leapfrog, which uses the above equations, would accelerate the path to the posterior \(\theta\). After coming up with proposal \(\theta\), we need to make a decision to reject or accept it. Therefore, \(\alpha\) would be created and compared to decide regarding accepting or rejecting the proposal. Here, instead of random walk for doing the proposal, we are solving for part of a particle, using the deterministic equation of motion. \(\alpha\) would be written as:  
\[\alpha=\frac{log(f(\theta_{2})-0.5m_{2}^{\intercal}M^{-1}m_{2}}{log(f(\theta_{1}) -0.5m_{1}^{\intercal}M^{-1}m_{1}} \tag{9}\]  
where subscript 2 is an updated version of parameters and 1 is the previous values. The above value in Equation (9) would be compared with \(u^{*}\), where \(u^{*}\) is a random uniform number between 0 and 1. The \(\alpha\) equation is intuitively highlighting that when \(f(\theta_{i})>f(\theta_{i-1})\), the proposal has a higher probability, and we accept the proposal. Therefore, if \(\alpha\) is greater than \(u^{*}\), we would accept the proposed value of \(\theta\) and m to be our next sample from our distribution, otherwise the value would remain as the previous values. In other words, if \(\alpha>u\) holds in Equation (9), then:  
\[m^{*}\leftarrow-m_{2}\ \&\ \theta^{*}\ \leftarrow\ \theta_{2}\ \ \text{ Otherwise }m^{*}\gets m_{1}\ \&\ \theta^{*}\leftarrow\ \theta_{1}. \tag{10}\]  
where \(m^{*}\) and \(\theta^{*}\) are proposal values, which would be used in the next iteration. We flip the sign of momentum in Equation 10 as if we start at a particular value of \(\theta\) and \(m\), then the part of deterministically leads to \(\theta^{*}\)and \(m^{*}\), or \(P(\theta^{*},m^{*}|\theta,m)=1\), but the reverse probability is 0. Therefore, what we do is switch the proposed \(m\) with its negative.