<\/p>\n\n\n\n
Probabilistic graphical models (Lauritzen, 1996) have garnered increasing interest in network psychometrics<\/a>, finding applications in diverse fields such as mental health (Hallquist et al., 2021), intelligence (Van Der Maas et al., 2017), and personality (Costantini et al., 2019). Many of these studies assume multivariate normality and employ Maximum Likelihood methods (Isvoranu & Epskamp, 2023) to fit Gaussian graphical models (GGMs). Some researchers have explored the use of Bayesian approaches (e.g., Williams & Mulder, 2020a) for GGM estimation. Despite the advantages of Bayesian models, some limitations of existing Bayesian GGM methods may hinder the development of the field. For instance, currently, different data types (e.g., binary or ordinal vs.\u00a0continuous data) require the development of new models and sampling schemes (e.g., Marsman et al., 2022). However, Bayesian GGMs can be formulated with more general and flexible model specifications. This study aims to introduce a general approach for regularized Bayesian Gaussian graphical models (GBGGMs) and to illustrate the performance of proposed models through a toy simulation and an empirical example.<\/p>\n\n\n\n A graph can be defined as a pair \\(G = (V,E)\\), where \\(V\\) represents a set of nodes and \\(E\\) the edges connecting such nodes. \\(E\\) is by itself a subset of the set \\(V \\times V\\)\u00a0which consists of pairs of distinct nodes, resulting in \\(G\\) having no cases with multiple edges or loops (Lauritzen, 1996). In a Gaussian Graphical Model (GGM), each node of the graph represents a normally distributed variable, and the set of observed variables in the model follow a multivariate normal distribution. Each edge in a GGM represents the partial correlation between each node, as the variables are conditionally dependent given the remaining variables in the network (Epskamp, Waldorp, et al., 2018). These partial correlations are obtained from the concentration matrix \\(\\Sigma^{-1}\\), which is the inverse of the covariance matrix \\(\\Sigma\\). From the elements \\(p\\) of the concentration matrix, the partial correlation for a pair of variables \\(i\\) and \\(j\\), \\(\\omega_{ij}\\), is defined as:<\/p>\n\n\n\n <\/p>\n\n\n\\( Compared to Pearson correlation, partial correlations between two variables have the advantage of being adjusted for the influence of other variables, thus avoiding spurious correlations. Additionally, under specific circumstances, partial correlations can indicate causal relationships between variables (Andersson et al., 1997; Epskamp, Maris, et al., 2018). In the graphical representation of a GGM, the absence of an edge between nodes indicates the absence of a relationship, whether marginal or conditional, between a pair of variables.<\/p>\n\n\n\n The application of network methods in psychological research has gained momentum over the past five years. According to Isvoranu and Epskamp (2023), several factors contribute to this trend, including the need for more intuitive and visual representations of psychological constructs and their interrelations, the exploration of potential causal relationships<\/a> in correlational studies, and efforts to address the limitations of reflective theories of mental constructs. Although the use of network methods in psychology has recently surged, its roots trace back to the 1930s when Moreno (1934) first proposed their use in studying social networks or sociograms. Moreno\u2019s pioneering work primarily focused on the social connections individuals establish. Consequently, networks predominantly depicted relationships between individuals rather than observed variables (Wasserman & Faust, 1994). In contrast, contemporary network psychometrics (Epskamp, Maris, et al., 2018) employs graphical models where nodes represent variables and edges represent statistical relationships between them (Isvoranu & Epskamp, 2023). This approach bears similarities to Structural Equation Models (SEM) but from a more exploratory perspective (Epskamp, 2020).<\/p>\n\n\n\n A crucial aspect of Gaussian Graphical Models (GGMs) involves identifying partial correlations that can be reliably inferred as exactly zero (Williams, 2020). Many estimated relationships in GGMs may be spurious correlations<\/a>, warranting a method to address this issue. Regularization methods, such as the graphical LASSO (or \u201cglasso\u201d), offer a solution (Friedman et al., 2008). The glasso technique involves applying a modified LASSO regression to each variable, which effectively forces small partial correlation coefficients to zero. Similar to LASSO regression, the glasso method utilizes a tuning parameter to control the degree of sparsity in the resulting graph. Lower values of this parameter yield denser graphs, while higher values produce sparser ones. Additionally, glasso can be combined with the Extended Bayesian Information Criterion (EBIC) fit statistic for model selection. This combined approach has proven to be effective in determining an optimal value for the tuning parameter (Isvoranu & Epskamp, 2023).<\/p>\n\n\n\n Despite the widespread use of glasso, Williams (2020) has argued that it: (i) leads to overshrinkage when the sample size is smaller than the number of variables; (ii) leads to overshrinkage when a set of predictors is correlated; (iii) has a higher prediction error than the ridge when the predictors are highly correlated; (iv) can lead to overshrinkage of large coefficients; and (v) does not perform well if the true data generating process is not known. Due to these limitations, alternative regularization methods for GGMs have been proposed by various authors (see Williams (2020) for a survey of alternative Maximum Likelihood methods). This scenario highlights the necessity of developing new regularization methods for GGMs, given that a specific method may fail in specific contexts. Despite of growing interest in Bayesian approaches form GGMs<\/a>, they still lag behind in terms of diversity compared to Maximum Likelihood-based methods.<\/p>\n\n\n\n In response to the limitations of regularization methods within the Maximum Likelihood approach, several authors have proposed Bayesian alternatives (e.g., Byrd et al., 2021; Li et al., 2019; Liu & Martin, 2019; Marlin & Murphy, 2009; Marsman & Haslbeck, 2023; A. Mohammadi & Wit, 2015; R. Mohammadi et al., 2021; Park et al., 2022; Sagar et al., 2024; Wang, 2012). In the context of network psychometrics, probably the most well-known approach is one proposed by Williams (2021). The author proposes that estimation follows a conventional Bayesian approach used for modeling multivariate normal variables by sampling the covariance matrix from an inverse Wishart prior (Sun & Berger, 2007). The Wishart distribution involves two hyperparameters: the base covariance matrix \\(V\\) and the degree of freedom \\(k\\), where \\(k\\)\u00a0is equal to or greater than the number of variables. This estimation approach aligns with the estimation of regularized Bayesian Gaussian graphical models<\/a> (BGGMs), wherein covariance or correlation values are shrunk, or regularized, toward zero. This is achieved by setting \\(V\\) as a diagonal (or identity) matrix, with \\(k\\) determining the degree of regularization in the model; higher values indicate stronger regularization. To discern the presence or absence of edges, post-hoc statistics based on the posterior distribution of edges are employed. A straightforward strategy entails examining whether zero is encompassed within the posterior distribution of estimated partial correlations. If so, the corresponding edge should be omitted from the final estimate of the GGM.<\/p>\n\n\n\n While the approach proposed by Williams (2021) has notable strengths, we argue that a generalized Bayesian GGM approach should possess certain key capabilities: (i) the ability to model correlations across diverse data types; (ii) the provision of parameter estimates indicating the probability of edge inclusion; and (iii) flexibility and ease of implementation with general-purpose optimization algorithms. For the best of our knowledge, these desiderata are not implemented concomitantly in any current Bayesian GGM method. In regards to William\u2019s approach, from this standpoint, we identify two primary limitations. Firstly, it is tailored specifically for continuous data. Although similar methods exist that employ extensions of the Gibbs sampling scheme (e.g., Albert & Chib, 1993; Hoff, 2007; Talhouk et al., 2012), and are incorporated into the BGGM R package (Williams & Mulder, 2020b), these approaches are heavily reliant on their particular sampling schemes, limiting their adaptability to alternative estimation methods such as different Markov Chain Monte Carlo (MCMC) algorithms or variational Bayes and Laplace Approximation. Secondly, Williams\u2019s approach does not directly estimate the presence or absence of edges during the estimation process. Instead, this information is derived post-estimation from the posterior distribution of parameters using procedures like the Savage-Dickey ratio (Williams & Mulder, 2020a).<\/p>\n\n\n\n One main challenge with the Bayesian estimation of GGMs is the sampling of positive semi-definite matrices. In the case of Markov Chain Monte Carlo (MCMC) estimation (Hanada & Matsuura, 2022), a sampler based on the Wishart distribution can be developed, with the base covariance matrix \\(V\\) set as the current sample of the covariance matrix, and the degrees of freedom \\(k\\)\u00a0representing the \u201cprecision\u201d of the sampler: larger values of \\(k\\)\u00a0will make the sampler draw covariance matrices closer to \\(V\\). One limitation of this approach is that it involves sampling the entire matrix simultaneously, which limit the precision of the posterior distributions of the estimated covariances to be equal. Also, the sampling of the entire matrix at once does not allow to fix a subset of parameters at known (or assumed to be known) values.<\/p>\n\n\n\n Additionally, in the scenario of point estimation using deterministic optimization methods like gradient descent, another challenge arises: all (co)variances must be simultaneously constrained to ensure the resulting estimated matrix is positive semi-definite. The Generalized Bayesian Gaussian Graphical Model (GBGGM) approach proposed in this paper, which is suitable both for MCMC as well as deterministic optimization, has as its initial step the use of the Cholesky decomposition approach for the estimation of the correlation parameters (Archakov et al., 2024; Archakov & Hansen, 2021; Bhat & Mondal, 2021).<\/p>\n\n\n\n The utilization of Cholesky decomposition for estimating correlation matrices is well-documented in the literature (the interested reader is referred to Chen & Leng, 2015; Ghosh et al., 2021; Roverato, 2000; Verzelen, 2010; Ye et al., 2020). In this section, we outline the fundamental steps for its Bayesian estimation. The Cholesky decomposition of matrix \\(A\\) is defined as:<\/p>\n\n\n\\( where \\(L\\) is a lower triangular matrix and \\(L^T\\)\u00a0is its transpose. This decomposition guarantees that \\(A\\) will be a positive semi-definite matrix, as long as the diagonal elements of \\(L\\) are equal or larger than 0. For the estimation of a correlation matrix, the lower-diagonal elements of \\(L\\) are treated separately and can be shifted (in the case of deterministic optimization), or sampled (in the case of the MCMC estimation), independently. However, to guarantee that \\(A\\) will be a correlation matrix \\(R\\), some additional steps are necessary. More specifically, starting with a lower triangular matrix \\(C\\), where all the diagonal elements are equal to 1 and the lower-diagonal elements, or the \u201cweights\u201d, are real, \\(L\\) is calculated with:<\/p>\n\n\n\\( where \\(E\\) is a diagonal matrix with its diagonal values equal to \\(e_i\\), for \\(i=1,2,3,\u2026,k\\)\u00a0is the order of \\(E\\), and \\(e_i\\) is equal to the inverse of the Euclidean norm of the i<\/em>-th row of \\(C\\). The estimation of correlation matrices with the Cholesky decomposition, then, can be conducted on the elements of \\(C\\), also called in this paper as \u201cweights\u201d, which are used to estimate the correlation and partial correlation matrices (a similar approach is also found in Lewandowski et al., 2009). To estimate a covariance matrix, this approach can be extended by using the LDL decomposition (Krishnamoorthy & Menon, 2013), which is defined as:<\/p>\n\n\n\\( where \\(D\\) is a positive diagonal matrix and \\(L\\) is the same as before.<\/p>\n\n\n\n In our generalized approach to BGGMs, the modeling and estimation process includes an additional step. If we define \\(Q\\) as the inverse of the lower-diagonal matrix \\(L\\), \\(Q=L^{-1}\\), and considering that the partial correlation matrix is derived from the inverse of the correlation matrix, we can exploit the fact (Krishnamoorthy & Menon, 2013) that<\/p>\n\n\n\\( to set a prior to the lower-diagonal elements of \\(Q\\), which will correspond to the weights of the partial correlation matrix. Then, to estimate a regularized partial correlation matrix, one can use shrinkage priors on the \\(p=k(k-1)\/2\\)\u00a0lower-diagonal elements (i.e., \\(q_j\\)) of \\(Q\\), as similarly performed in Bayesian penalized regression (van Erp, 2020). In summary, we set the same symmetric distribution as the prior for all the lower-diagonal elements of \\(Q\\), where this prior has to have its location parameter equal to 0 and its scale parameter serving as the inverse of the strength of the regularization (i.e., amount of shrinkage). It is also possible to use distributions with additional parameters, as long as there are at least a location and a scaling parameter.<\/p>\n\n\n\n From this, the general representation of the approach is the following:<\/p>\n\n\n\\( In this case, we are using \\(\\boldsymbol{\\Theta}\\) to represent the set of additional parameters (i.e., other possible parameters that are not a location parameter, which is always set to 0) that a distribution may have, including, at least, a scaling parameter. Therefore, for instance, ridge regularization (van Erp et al., 2019) is achieved by setting the prior \\(P(0,\\boldsymbol{\\Theta})\\)\u00a0to be the Normal distribution with a mean of 0 and standard deviation of \\(\\lambda\\). More specifically, we have<\/p>\n\n\n\\( where \\(\\lambda\\)\u00a0is the shrinkage parameter, \\(p\\) is the number of parameters, and \\(Cauchy^+\\)\u00a0is the half-Cauchy distribution. Bayesian LASSO regularization (van Erp et al., 2019), in contrast, is achieved with<\/p>\n\n\n\\( where Laplace<\/em> stands for the Laplace, or double-exponential, distribution.<\/p>\n\n\n\n There are several generalizations and extensions of the ridge and LASSO methods in the Bayesian literature (van Erp et al., 2019). In general, any symmetric distribution centered at zero (D\u00e9niz, 2021) can be used as the prior for the elements of \\(Q\\). For instance, one could use the logistic or the Cauchy distributions instead of the normal or the Laplace distributions. One could also add an extra parameter to the model and then set the prior to be equal to the three-parameter \\(t\\) distribution (Lange et al., 1989), a heavy-tailed version of the ridge, or to the double lomax distribution (Bindu & Sangita, 2015), a heavy-tailed version of the LASSO. Further still, it is possible to set the prior as less \u201cwell-behaved\u201d distributions, such as the generalized normal (Nadarajah, 2005), the skew-normal (Arellano-Valle & Azzalini, 2006), or the Normal-exponential-gamma (Scheipl & Kneib, 2009) distributions. However, caution should be taken when choosing these priors, as the shrinkage properties of these alternatives may not be appropriate for the objective at hand. As argued before, in regards to the LASSO regularization, the performance of different regularization methods can be sub-optimal.<\/p>\n\n\n\n After setting the prior for the elements of \\(Q\\), the next step is to define the likelihood function that represents the data. For continuous data, the likelihood function is typically set as the multivariate normal distribution. For example, when dealing exclusively with continuous data and aiming to implement a correlation model with Bayesian LASSO regularization, the complete model can be specified as follows<\/p>\n\n\n\n <\/p>\n\n\n\\( In psychometrics, a critical modeling challenge involves estimating polychoric correlations (J\u00f6reskog, 1994). Polychoric correlations posit that the observed monotonic relationships among ordered categorical variables stem from discretizing linearly related Gaussian variables with unobserved thresholds. This assumption permits estimation of joint response proportions for ordered categorical variables based on the probability function of a multivariate normal distribution. However, the full Maximum Likelihood method is often impractical due to the exponential growth in joint responses with the number of variables and response categories (Olsson, 1979). An alternative feasible method within the Maximum Likelihood framework is to estimate polychoric correlations from bivariate normal marginal distributions given thresholds (J\u00f6reskog, 1994). Yet, this method often yields a correlation matrix that is not positive definite. In contrast, our proposed generalized approach ensures the estimated correlation matrix is positive definite while benefiting from the computational efficiency of calculating likelihood from bivariate normal marginal distributions. Additionally, our method enables simultaneous estimation of thresholds without necessitating a Bayesian two-step alternative (Lee & Poon, 1987).<\/p>\n\n\n\n Formally, the logarithm of the likelihood for polychoric correlations between each pair of variables \\(i\\) and \\(j\\) is defined as \\(\\text{ln}L_{ij}=N\\sum_{a=1}^{M_{ij}} p_{ija}\\text{ln}\\pi_{ija}(r_{ij}, \\tau_{ija})\\). Here, <\/em>\\(N\\) is the sample size and \\(M_{ij}\\)\u00a0represents the total number of possible combinations of response categories for variables \\(i\\) and \\(j\\). For instance, if variable \\(i\\) has 3 response categories and \\(j\\) has 4\u00a0response categories, then \\(M_{ij} = 12\\). \\(p_{ija}\\)\u00a0denotes the proportion of responses in the combination \\(a\\) of response categories for variables \\(i\\) and \\(j\\), and \\(pi_{ija}(r_{ij}, \\tau_{ija})\\)\u00a0is the cumulative probability of the bivariate normal distribution with correlation \\(r_{ij}\\)\u00a0and threshold parameter \\(\\tau_{ija}\\). For simplicity, we define \\(PolyBVN(\\hat{R}, \\boldsymbol{\\tau})\\)\u00a0as the polychoric likelihood based on the bivariate normal distribution, where \\(\\hat{R}\\) is the estimated correlation matrix and \\(\\boldsymbol{\\tau}\\)\u00a0is the vector of threshold parameters for all the items. The final model can be specified as follows:<\/p>\n\n\n\\( As a key feature of our proposed generalized Bayesian Gaussian Graphical Model method, we aimed to incorporate the estimation of network sparsity as part of the model per se. In previous Bayesian GGM methods found in the literature, sparsity of the correlation matrix was typically estimated through post-hoc procedures conducted after model fitting. In contrast, we propose utilizing a mixture probability model akin to the approach suggested by Gorbach et al. (2020). The main idea with this approach is that the off-diagonal elements of the partial correlation matrix implied by the model form groups of \u2018connected\u2019 and \u2018disconnected\u2019 nodes. To achieve this, the extension of the model requires the calculation of the implied partial correlation matrix \\(\\hat{\\Omega}\\)\u00a0and some additional parameters.<\/p>\n\n\n\n While complex models like the one proposed by Gorbach et al. (2020) may enhance estimation quality, we advocate for a simpler approach. We exemplify our proposition with the same example as before: imagine we are dealing with continuous data and seeking to implement a correlation model with Bayesian LASSO regularization. If one desires to estimate a sparse partial correlation matrix, the complete model can be specified as follows<\/p>\n\n\n\\( The initial step in the extension proposed within this model involves computing \\(\\delta_j\\), representing the probability of the edge \\(j\\)\u00a0being present. This probability is calculated as a scaled transformation of the partial correlation \\(\\omega_j\\). Specifically, \\(F(\\cdot)\\)\u00a0is the Fisher transformation (i.e., the arc-tangent function) and \\(\\theta\\)\u00a0is a scaling parameter. \\(H\\phi\\)\u00a0is the cumulative function of the half-normal distribution, which depends on the parameter \\(\\sigma\\). Under this setup, larger values of \\(\\theta\\)\u00a0increase the slope of the relation between the magnitude of the partial correlation and the probability of presence of the edge. Larger values of \\(\\sigma\\) diminish the maximum possible probability of edge presence across all edges. Subsequently, the likelihood (or, probably more adequately, a hyperprior) of the Fisher-transformed partial correlations \\(F(\\omega_j)\\)\u00a0is calculated using a mixture distribution. The mixture comprises a normal density, with mean 0\u00a0and standard deviation defined by the \\(\\gamma\\)\u00a0parameter, weighted by \\((1-\\delta_j)\\), and a Bernoulli mass, with the probability parameter equal to \\(\\delta_j\\), weighted by \\(\\delta_j\\). The scaling parameters \\(\\theta\\), \\(\\sigma\\), and \\(\\gamma\\) are assigned half-Cauchy priors.<\/p>\n\n\n\n This extended version of the basic generalized Bayesian Gaussian graphical model introduces an alternative method for determining whether a node should be included in the final estimation of the partial correlation matrix. Notably, this approach eliminates the need for post-estimation methods to deduce the presence or absence of an edge from the posterior distribution of parameters. Instead, this model enables the use of point estimation, rather than full posterior distribution estimation, to make conclusions regarding edge presence. While the richness of Bayesian inference typically arises from the posterior distribution, in certain applications, the Maximum A Posteriori (MAP) estimate may suffice. Also, the post-estimation methods can still be used with the resulting estimates, if the full posterior distribution is estimated. Therefore, our model adds flexibility, but also additional information that can be used in the decision process about the true network that generates the data.<\/p>\n\n\n\n In summary, our suggestion for a generalized approach in estimation of BGGMs consists of the following steps:<\/p>\n\n\n\n 1. Sample the elements of \\(C\\).<\/p>\n\n\n\n 2. With C<\/em>, estimate \\(L\\), \\(Q\\), \\(R\\), and \\(\\Omega\\).<\/p>\n\n\n\n 3. Choose a prior distribution for the lower-diagonal elements \\(q_j\\)\u00a0of \\(Q\\).<\/p>\n\n\n\n 4. Choose any extra necessary shrinkage parameters according to sparsity constraints.<\/p>\n\n\n\n 5. Calculate the likelihood of your data.<\/p>\n\n\n\n In this section, we conduct a simple toy simulation to provide an initial assessment of the feasibility of our approach. To achieve this, we simulated a network with 7 nodes (equating to 21 correlations) based on a random Directed Acyclic Graph [DAG; Pearl (2009)], using the \u201cordered\u201d method from the bnlearn R package (Scutari & Silander, 2024). From the random DAG, we derived its moral graph (i.e., the undirected graph that represents the equivalence class of the DAG) and, for the included edges, we generated random values using a random distribution with lower bound of .30 and upper bound of .70. If the values did not allow the transformation of this base matrix (i.e., the partial correlation matrix) to a correlation matrix, we repeated the sampling procedure until this transformation was feasible. The resulting partial correlation matrix is presented in Table 1.<\/p>\n\n\n\n Table 1<\/strong><\/p>\n\n\n\n The true simulated partial correlation matrix<\/em><\/p>\n\n\n\nRegularized Gaussian Graphical Models (GGMs)<\/h3>\n\n\n\n
\n\\begin{equation}
\n\\tag{1}
\n\\omega_{ij} = -\\frac{p_{ij}}{\\sqrt{p_{ii}p_{jj}}}
\n\\end{equation}
\n\\)\n\n\n\nThe Standard Bayesian Estimation of GGMs: Williams (2021)<\/h3>\n\n\n\n
Cholesky Decomposition of Correlation Matrices<\/h3>\n\n\n\n
\n\\begin{equation}
\n\\tag{2}
\nA = LL^{T}
\n\\end{equation}
\n\\)\n\n\n\n
\n\\begin{equation}
\n\\tag{3}
\nL = (C^{T}E)^{T}
\n\\end{equation}
\n\\)\n\n\n\n
\n\\begin{equation}
\n\\tag{4}
\nA = LDL^{T}
\n\\end{equation}
\n\\)\n\n\n\nA General Approach for Bayesian Estimation of GGMs<\/h3>\n\n\n\n
\n\\begin{equation}
\n\\tag{5}
\nA^{-1} = Q^{T}Q
\n\\end{equation}
\n\\)\n\n\n\n
\n\\begin{gather*}
\n\\tag{6}
\nq_j \\sim P\\left(0,\\boldsymbol{\\Theta}\\right)\\text{, for } j=1,\u2026,p\\text{,}\\\\
\n\\end{gather*}
\n\\)\n\n\n\n
\n\\begin{gather*}
\n\\tag{7}
\nq_j \\sim N\\left(0,\\frac{1}{\\lambda}\\right)\\text{, for } j=1,\u2026,p\\text{,}\\\\
\n\\lambda \\sim Cauchy^{+}(0,1)
\n\\end{gather*}
\n\\)\n\n\n\n
\n\\begin{gather*}
\n\\tag{8}
\nq_j \\sim Laplace\\left(0,\\frac{1}{\\lambda}\\right)\\text{, for } j=1,\u2026,p\\text{,}\\\\
\n\\lambda \\sim Cauchy^{+}(0,1)
\n\\end{gather*}
\n\\)\n\n\n\nDefining the Likelihood of the GBGGM<\/h3>\n\n\n\n
\n\\begin{gather*}
\n\\tag{9}
\nq_j \\sim Laplace\\left(0,\\frac{1}{\\lambda}\\right)\\text{, for } j=1,\u2026,p\\text{,}\\\\
\n\\lambda \\sim Cauchy^{+}(0,1)\\\\
\n\\hat{L} \\leftarrow (\\hat{C}^{T}\\hat{E})^{T}\\\\
\n\\hat{R} \\leftarrow \\hat{L}\\hat{L}^{T}\\\\
\nX \\sim MVN(\\boldsymbol{0}, \\hat{R})
\n\\end{gather*}
\n\\)\n\n\n\n
\n\\begin{gather*}
\n\\tag{10}
\nq_j \\sim Laplace\\left(0,\\frac{1}{\\lambda}\\right)\\text{, for } j=1,\u2026,p\\text{,}\\\\
\n\\lambda \\sim Cauchy^{+}(0,1)\\\\
\n\\hat{L} \\leftarrow (\\hat{C}^{T}\\hat{E})^{T}\\\\
\n\\hat{R} \\leftarrow \\hat{L}\\hat{L}^{T}\\\\
\nX \\sim PolyBVN(\\hat{R}, \\boldsymbol{\\tau})\\\\
\n\\boldsymbol{\\tau} \\sim N(0,1)
\n\\end{gather*}
\n\\)\n\n\n\nEstimating the Sparsity of the Network<\/h3>\n\n\n\n
\n\\begin{gather*}
\n\\tag{11}
\nq_j \\sim Laplace\\left(0,\\frac{1}{\\lambda}\\right)\\text{, for } j=1,\u2026,p\\text{,}\\\\
\n\\lambda \\sim Cauchy^{+}(0,1)\\\\
\n\\hat{L} \\leftarrow (\\hat{C}^{T}\\hat{E})^{T}\\\\
\n\\hat{R} \\leftarrow \\hat{L}\\hat{L}^{T}\\\\
\nX \\sim MVN(\\boldsymbol{0}, \\hat{R})\\\\
\n\\\\
\n\\delta_j \\leftarrow H\\phi(|\\theta \\cdot F(\\omega_j)|, \\sigma)\\text{, for } j=1,\u2026,p\\text{,}\\\\
\n\\theta \\sim Cauchy^{+}(0,1)\\\\
\n\\sigma \\sim Cauchy^{+}(0,1)\\\\
\nF(\\omega_j) \\sim (1-\\delta_j)N(0, \\gamma) + \\delta_{j}Bernoulli(\\delta_j)\\text{, for } j=1,\u2026,p\\text{,}\\\\
\n\\gamma \\sim Cauchy^{+}(0,1)
\n\\end{gather*}
\n\\)\n\n\n\nIllustration With a Toy Simulation<\/h2>\n\n\n\n