## 4 Empirical Study
### Results with LGPIF (H = 1679)  
We construct a training set with 1276 records for modeling and a test set with 403 records for validation. As shown in Figure 4, Model (A) serves as the gold standard, Model (B) represents a flawed model that reflects real-world risks, Model (C) applies the Gustafson correction, and Model (D) utilizes the SIMEX correction within a traditional GLM framework. We systematically evaluate the performance of Model (B), -(C), and -(D) across three scenarios with varying error rates \(R_{\varepsilon_{x}}=0.01,R_{\varepsilon_{x}}=0.1\), and \(R_{\varepsilon_{x}}=0.4\) in the NDB covariate \(\mathbf{x}^{5*}\). Additionally, we compare the effects of the scaling factor \(0<\zeta<1\) in Model (C) and -(D) within each error rate scenario. For each hierarchical GLM model, Model (A), -(B), and -(C) shown in Figure 4, we run two independent Markov chains with \(M=\) 60,000 iterations of Gibbs sampling, following Algorithm S2 in Part 2 of the Supplementary Materials. Since Model (A) is the gold standard, the focus is on comparing Model (B), -(C), and -(D) across datasets with error rates \(R_{\varepsilon_{x}}=0.01,R_{\varepsilon_{x}}=0.1\), and \(R_{\varepsilon_{x}}=0.4\). Model (C) employs the Gustafson correction with scaling factors \(\zeta=\{0.1,\ldots,0.9\}\). Within each Gibbs sampling iteration in Model (A), -(B), and -(C), a Metropolis-Hastings (MH) technique is embedded to update the outcome parameters--\(\beta_{j}^{F},\psi_{j},\beta_{j}^{S},\sigma_{j}^{2}\)--as conjugate priors are not available. The initial 10,000 iterations are treated as burn-in and discarded. Convergence is confirmed using the Brooks-Gelman statistic [36], ensuring the adequate mixing of the chains. For Model (D), a conventional GLM, parameters are estimated using the Maximum Likelihood Estimation (MLE) method. The aggregate claim amount model \(f(S_{h}|\mathbf{X}^{F},\mathbf{X}^{S},\xi_{h},\psi,\mu_{h},\sigma^{2})\), representing the risk premium function, is derived by integrating these two gold standard models through Monte Carlo simulation: the (i) estimated claim count model \(f(N_{h}|\mathbf{X}^{F},\xi_{h},\psi)\) and the (ii) estimated claim amount model \(f(\tilde{Y}_{h}|\mathbf{X}^{S},\mu_{h},\sigma^{2})\) (detailed results for each gold standard model are provided in Part 4 of the Supplementary Materials). The integration process begins by simulating a claim count sample \(N_{h}\) for a policy \(h\). This sample is drawn from the claim count model \(f(N_{h}|\mathbf{X}^{F},\xi_{h},\psi)\). Subsequently, for each drawn claim count \(N_{h}\), a corresponding set of claim amount samples is generated. These samples are obtained from the claim amount model \(f(\tilde{Y}_{h}|\mathbf{X}^{S},\mu_{h},\sigma^{2})\). This iterative process is repeated numerous times, allowing for the construction of a comprehensive distribution of aggregate claims \(S_{h}\) for a policy \(h\)[37]. The resulting predictive distribution of the aggregate claim amount modeling is illustrated in Figure 5.  
This figure compares the cluster-wise distributions of the aggregate claim amount \(S_{h}\) on a log scale with the overall distribution of \(S_{h}\). This comparison facilitates the evaluation of each cluster's risk profile and its impact on the total aggregate claim. Significantly, the distribution for 'Cluster 2' (shown as a dotted blue curve) closely resembles the shape of the overall distribution, indicating a substantial contribution to the overall risk premium. In what follows, we will examine the results of Model (A), -(B), -(C), and -(D), focusing exclusively on the claim amount component \(\tilde{Y}_{h}|\mathbf{X}^{S*}\sim\mathbf{LogN}(\mu_{h}=ln\big{(}E[\tilde{Y} _{h}|\mathbf{X}^{S*}]\big{)}-\frac{1}{2}\sigma_{j}^{2},\ \sigma_{j}^{2})\) (log-normal GLM).