## 4 Numerical Results  
In this section, we validate the model against numerical simulations, and we present some numerical results of the variation of the 90 th percentile of the exposure against different network parameters. In Figure 5, we validate the CDF from with \(\varphi(s)\) given by Theorem 1, with Monte-Carlo simulations of the same scenario, to verify the model we developed while the parameters used in the validation are presented in Table 4.  
We also perform a comparison between the newly created model and the model in our previous work [25] versus the exposure simulation derived directly from the NYUSIM data in Figure 8. This comparison shows that the old model overestimates the exposure especially at the lower percentiles. This overestimation can be attributed to the fact that the old model assumes uniform array pattern, which is equivalent of having an isotropic antenna gain at the transmitter, and gain values towards the MT. It is also apparent, in the model derived in this study compared to the simulations, the small error between the simulation and our model which can be attributed to the error in fitting finite gain data into infinite distributions alongside the errors from the imperfect fitting. The former error can be calculated depending on the maximum gain the antenna can produce as in Equation (16) where \(f_{X}(x)\) is the probability density function of the fitted distribution, and \(G_{max}\) is the maximum gain the antenna array can produce. \[er_{fit}=\int_{G_{max}}^{\infty}f_{X}(x)dx \tag{16}\]  
Since massive MIMO will be implemented in really diverse scenarios, it is important to investigate the effect of each variable on the exposure. For this, we perform a variance-based sensitivity analysis on the calculated model which estimates the effect of the variance of the inputs on the variance of the output as form of indices called Sobol indices [32]. Sensitivity analysis is often performed to determine the importance of the input variables on the output of the model. \(n\)th degree Sobol indices \(S_{1,2...n}\) determine the fraction of the output variance attributed to the set of inputs of degree \(n\). We are interested here in determining the total sobol indices which can be defined as in Equation (17). The total Sobol indices determine the variance of the output due to the variance of an input, in addition to the \begin{table}
\begin{tabular}{c c} \hline \hline
**Parameter** & **Value** \\ \hline \(b\) & 1.45 \\ \(c\) & 966.5 \\ \(a\) & 0.57 \\ \(\lambda\) & \(2\times 10^{-5}\) \\ \(\alpha\) & 0.5 \\ \(\eta\) & 3 \\ \hline \hline \end{tabular}
\end{table} 
Table 5: Parameters used in the analysis simulations for Figures 6 and 7. Figure 7: 90 th percentile of the exposure as function of system utilization for different values of the BS density. interactions between the specified input variable and the other input variables. It is an alternative to computing the higher-order sobol indices for every variable. \[S_{Ti}=\frac{E_{X_{\_},\_}\big{(}Var_{\_}X_{\_}i}(Y|X_{\_}{\_}{-i})\big{)}{Var(Y)} \tag{17}\]  
Since obtaining the exposure from the CDF requires solving an inverse function, it takes a long time to determine the sobol indices and it may introduce inaccuracies at extreme input values. To avoid this, we use polynomial chaos expantion (PCE) using Latin Hypercube samples to estimate a metamodel to represent the 90 th percentile of the exposure in the cell. PCE approximates the relation between the model's output to its inputs by expanding it in an orthogonal polynomial basis [32]. The metamodel can be denoted as  
\[\Upsilon(\zeta)=\sum_{\_}{\beta\in\mathbb{N}^{d}}a_{\_}\Psi_{\_}\beta(\zeta), \tag{18}\]  
where \(\Psi_{\_}\beta\) are multivariate polynomials defined as the product of univariate polynomials of degrees \((\beta_{\_}1,\ldots,\beta_{\_}d)\) and \(a_{\_}\beta\) are the polynomial's coefficients. The fit model has a leave-one-out error of 0.14%.