Computational framework for numerical probabilistic seismic hazard analysis (PSHA)

Abstract

A computer-implemented method of performing probabilistic seismic hazard analysis (PSHA) includes receiving a seismic source model, a ground motion model, and a probability distribution, constructing a joint probability distribution, using adaptive importance sampling to obtain an optimal sampling density, computing a PSHA hazard estimate and a hazard deaggregation distribution, and providing the PSHA hazard estimate and the hazard deaggregation to a user. Another computer-implemented method includes estimating a mean and fractile PSHA hazards for each of a series of iterations by drawing a random sample set from a proposal sampling density, evaluating an indicator function for each sample and multiply by an earthquake occurrence rate and a likelihood ratio to produce a weight vector, and estimating a mean hazard from an average of the weight vector.

Description

Patent Application U Cal No. BK-2024-147-2-PCT MN No.407869-0211 COMPUTATIONAL FRAMEWORK FOR NUMERICAL PROBABILISTIC SEISMIC HAZARD ANALYSIS (PSHA) TECHNICAL FIELD

[0001] This disclosure relates to probabilistic seismic hazard analysis (PSHA), more particularly to PSHA using adaptive intelligent sampling. BACKGROUND

[0002] Probabilistic Seismic Hazard Analysis (PSHA) has become a foundational method for determining seismic design levels and conducting regional seismic risk analyses since its first inception. The development of PSHA was driven by the necessity for a probabilistic framework to accurately quantify seismic hazards, acknowledging the unpredictable nature of earthquake location, magnitudes, and ground motion intensities. In PSHA, earthquake location, magnitude, and ground motions are treated as random variables, facilitating the computation of annual exceedance probabilities at various ground motion intensities. Integrating hazard with fragility curves further enables the determination of annual probabilities of structural failures, thereby underscoring the methodology’s critical role in risk assessment.

[0003] Past data of earthquakes that include date, time, location, magnitude, ground motion, among many others, has been compiled into databases referred to as earthquake catalogs that include massive amounts of data. One can generate these catalogs to provide new data sets with particular characteristics. Using this data to make accurate seismic hazard predictions comprises a multi-faceted problem because of the complexity and the amount of the data. Since one cannot solve PSHA analytically due to the complexity in seismic source and ground motion models, numerous researchers have developed computer software for PSHA. Recently, the PSHA computer code verification project furnished insights into other software’s accuracy and provided benchmark hazard curves for various numerical examples. Existing software primarily employs Riemann summation for numerical integration of PSHA, which partitions the earthquake magnitude, location, and ground motion random variable into fine grids to approximate the actual integration. The Riemann summation offers robust PSHA integration with sufficiently dense grids. However, this method generally incurs a significantcomputational load in multi-dimensional integrations, exponentially increasing with the number of grids and dimensions. Furthermore, the results are highly sensitive to the chosen grid design (e.g., the initial point of the grid, grid spacing), leading to significant deviations from one code to another, especially for low exceedance probabilities.

[0004] Alternatively, other software adopted Monte-Carlo (MC) integration for PSHA. MC integration calculates exceedance probabilities by generating synthetic earthquake catalogs based on the seismic source and ground motion models and evaluating the recurrence of various ground motion intensities. MC integration’s primary advantage lies in its straightforward concept compared to Riemann summation, without the need to divide the integration range into small slices. Nevertheless, under MC framework, accurately estimating hazard from rare events requires a substantially long synthetic catalog, making the computation costly, especially for large ground motions with low exceedance probabilities, e.g., ^^^^ < 10−4 / yr.

[0005] Importance Sampling (IS) can offer a solution to the rare event simulation. IS was initially introduced in statistical physics to improve the computational efficiency of rare event simulation that would otherwise require a large sample size with conventional MC. IS relies on identifying an appropriate probability distribution (“IS distribution”) to explore low- probability spaces effectively. However, finding such an appropriate distribution can be challenging because there is no optimal sampling density that is universally applicable; rather, the proper selection of sampling density depends on the problem being solved. Thus, expensive numerical experiments are often conducted first through trial-and-error to identify IS distributions.

[0006] In regional seismic risk analysis, numerous studies have been conducted to sample hazard-consistent earthquake ground motions, some approaches first reduced the computational burden of seismic hazard and risk by defining IS distributions that samples large-magnitude earthquakes with a high probability. One approach expanded the approach by defining IS distributions to sample high-intensity ground motions. However, this approach highlighted the computational challenges to identify an effective IS distribution and ended up using K-mean clustering to reduce the number of ground motion samples. Another approach applied system reliability methods to calculate PSHA, and they selected the IS sampling density as the normal distribution centered at the "design point" derived from the first- and second-order reliability method.

[0007] Current approaches are computing intensive, especially when calculating lower probability PSHA. BRIEF DESCRIPTION OF THE DRAWINGS

[0008] FIG.1 shows a graphical representation of non-parametric importance sample (VEGAS) iterations of intelligent sampling (IS) in the adaptive importance sampling (AIS) framework of the embodiments.

[0009] FIG.2 shows seismic source geometry for three examples used in testing the AIS framework of the embodiments.

[0010] FIG.3 shows a graph of benchmark Probabilistic Seismic Hazard Analysis (PSHA) form the examples of FIG.2

[0011] FIG.4 shows graphs of comparisons of computational time for PSHA for different methods and different ground motions.

[0012] FIG.5 shows standard deviations of Monte Carlo (MC) estimates for conventional MC and AIS.

[0013] FIGs.6A-6D shows box plots showing the distribution of AIS MC estimates with different numbers of samples.

[0014] FIG.7 shows a comparison of the convergence of earthquake magnitude, source-to- distance, and the standard normal random variable for generating earthquake ground motion between IS densities and marginal distributions of hazard deaggregation at different ground motion intensities using 10,000 samples.

[0015] FIG.8 shows Kolmogorov-Smirnov D statistic and mean differences between hazard deaggregation and the proposed optimal sampling density for an areal source.

[0016] FIG.9 shows standard deviations of the conventional MC, IS MC and AIS estimates as a function of computation time at different ground motions.

[0017] FIGs.10A-10D shows a box plot showing the distribution of AIS MC estimates for a fault source with different numbers of samples.

[0018] FIG.11 convergence of earthquake magnitude, source-to-distance, and the standard normal random variable for generating earthquake ground motion between IS densities and marginal distributions of hazard deaggregation at different ground motion intensities using 1,000,000 samples.

[0019] FIG.12 shows Kolmogorov-Smirnov D statistics and mean differences between hazard deaggregation and the proposed optimal sampling density for a fault source.

[0020] FIG.13 shows a comparison of the standard deviations for convention MC, IS MC, and full AIS MC and partial AIS MC for combined sources.

[0021] FIGs.14A-14D shows a box plot of the distribution of full AIS MC estimates with different numbers of samples.

[0022] FIGs.15A-15D shows a box plot of the distribution of partial AIS MC estimates with different numbers of samples.

[0023] FIG.16 shows a flow chart of a method performing PSHA using the AIS approach of the embodiments.

[0024] FIG.17A shows an example of geometry of a seismic source and site location, and FIG.17B shows a continuous distribution of epistemic uncertainty variables.

[0025] FIGs.18A-18C show mean, 16thand 84 fractile hazard curves, and FIGs.18D-18F shows relative errors in the hazard curves.

[0026] FIG.19 shows a diagram of operation of an embodiment of a PSHA framework.

[0027] FIGs.20A-20C shows a seismic source geometries for an area source and a fault source, and aleatory and epistemic uncertainty variables.

[0028] FIGs.21A-21B show coefficients of variation for Monte Carlo results and Gaussian Population Monte Carlo Adaptive Importance Sampling results.

[0029] FIGs.22A-22B show Kolmogorov-Smirnov D statistics of the hazard distribution as a function of a total number of estimated hazards for areal and fault source examples.

[0030] FIG.23 shows sensitivity analysis for an areal source case for the G-PMC-AIS embodiment. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0031] Probabilistic Seismic Hazard Analysis (PSHA) relies on two widely utilized approaches with high computational demands: (a) Riemann sum and (b) conventional Monte Carlo (MC) integration. The first requires sufficiently fine slices across magnitude, m, distance, and ground motion, X, and the second requires extensive synthetic earthquake catalogs to compute seismic hazards accurately. These approaches are notably resource-intensive for low-probability seismic hazards, e.g., up to 108MC samples for a hazard with 10−4probability to achieve COV (Coefficient of Variance, or Covariance) of 1%. The term “earthquake catalog” refers to a database or databases that has a compilation of data bout past earthquakes. An earthquake catalog includes data about the time, date, location, magnitude, ground motion, and other parameters. These catalogs contain massive amounts of complex data and issues arise in determining the best way to sample the data. One can generate other catalogs from the databases of earthquake information to test processes for identifying earthquake hazards.

[0032] In PSHA, earthquake location, magnitude, and ground motions are treated as random variables, facilitating the computation of annual exceedance probabilities at various ground motion intensities. As used here, the term “exceedance” means an event where a ground motion parameter exceeds a specified design value during an earthquake, essentially that the intensity of a quake at a specific point is stronger than what structures at the point were designed to withstand.

[0033] At a site of interest, the annual frequency of ground motion exceedance from a single source can be calculated as:where ^^^^(^^^^ > ^^^^) is annual rate that ground motion, ^^^^, exceeds the target ground motion intensity, ^^^^. ^^^^ is the annual rate of earthquake occurrence greater than ^^^^minfrom the source, ^^^^ is earthquake magnitude, ^^^^minand ^^^^maxare minimum and maximum magnitudes considered, ^^^^ is the source-to-site distance, ^^^^minand ^^^^maxare minimum and maximum source- site distances, ^^^^ is standard normal random variable for generating earthquake ground motion, ^^^^min and ^^^^max are minimum and maximum ^^^^ (generally, ^^^^max ≥ 6 and ^^^^min ≤ −6 ), ^^^^^^^^,^^^^,^^^^(^^^^, ^^^^, ^^^^) is joint probability density function (PDF) of ^^^^, ^^^^, and ^^^^, ^^^^(^^^^ > ^^^^|^^^^, ^^^^, ^^^^) is the indicator function that takes 1 when ^^^^ > ^^^^, otherwise, 0. The ground motion ^^^^ given ^^^^, ^^^^, ℰ is generally calculated using ground motion models. The models usually assume the log normal distribution for the ground motion given explanatory variables such as ^^^^ and ^^^^. Naturally, these models provide the mean and standard deviation of logarithmic ground motion. Thus, the random ground motion can be calculated as: log X = μ(M, R) + ℰσ(M, R) (2)where ^^^^ and ^^^^ are mean and standard deviation of logarithmic earthquake ground motion. By taking exponential on both sides of Eq. (2), the ground motion ^^^^ can be calculated as X= ^^^^μ(M,R)+ℰσ(M,R)

[0034] If one assumes that the ground motion random variable ^^^^ is independent with respect to the ^^^^ and ^^^^ (Eq. (1)) can be modified aswhere ^^^^^^^^(^^^^), ^^^^^^^^|^^^^(^^^^|^^^^), ^^^^ℰ(^^^^) are PDF of ^^^^, ^^^^, and ^^^^. Under point source assumption, the distance ^^^^ and magnitude ^^^^ become independent random variables. Thus, the seismic hazard is given byThe total seismic hazard from multiple seismic sources (e.g., different faults) is the sum of each. Thus,where Λ(^^^^ > ^^^^) is the total annual frequency of exceedance of ground motion, ^^^^, ^^^^ is index for seismic sources, and ^^^^ is the total number of seismic sources. Under the assumption of Poisson process, the annual probability of exceedance, ^^^^(^^^^ > ^^^^), can be converted from annual frequency of exceedance (Λ; Eq. (5)) as ^^^^(^^^^ > ^^^^) = 1 − ^^^^−Λ(^^^^>^^^^) (6)

[0035] The embodiments can also disaggregate the total hazard (Eq. (6)) to better understand the earthquakes that contribute most to the hazard. Deaggregation, or breaking a particular hazard into its components, is also used to develop the select seismic records, such as from the earthquakes that contribute most to the hazard and conduct non-linear time-history analyses for the design of many critical buildings (Bazzurro and Allin Cornell, 1999; U.S. Nuclear Regulatory Commission, 2007). Mathematically, deaggregation of the hazard is the joint probability distribution of the ^^^^, ^^^^, and ^^^^ conditional on different levels of hazards ^^^^ toquantify the contributions of each component. The deaggregation of PSHA can be formulated using Bayes’ theorem as ^^^^(^^^^ > ^^^^ ∩ ^^^^ )( ^^^| , ^^^^, ^^^^^^^^, ^^^^, ^ ^^^^ > ^^^^) =^^^^(^^^^ > ^^^^)where ^^^^(^^^^ > ^^^^|^^^^, ^^^^, ^^^^) is the probability of ground motion ^^^^ is greater than ^^^^ given ^^^^, ^^^^, and ^^^^, ^^^^(^^^^, ^^^^, ^^^^) is joint probability of ^^^^, ^^^^, and ^^^^, and ^^^^(^^^^ > ^^^^) is the total probability that the ground motion is greater than ^^^^, which is the summation of ^^^^(^^^^ > ^^^^|^^^^, ^^^^, ^^^^)^^^^(^^^^, ^^^^, ^^^^) over all ^^^^, ^^^^, and ^^^^. Note that ^^^^(^^^^ > ^^^^|^^^^, ^^^^, ^^^^) can be expressed as an indicator function, ^^^^(^^^^ > ^^^^|^^^^, ^^^^, ^^^^), because the probability of ground motion ^^^^ greater than ^^^^ can only be 1 or 0 given ^^^^, ^^^^, and ^^^^ (see Eq. (2)).

[0036] By replacing probability mass function, ^^^^(^^^^, ^^^^, ^^^^), with probability density function, ^^^^^^^^,^^^^,ℰ(^^^^, ^^^^, ^^^^), and change the summation into integration, Eq. (7) can be expressed as:

[0037] By Eq. (1), the denominator of Eq. (8) equals ^^^^∕^^^^. Thus,

[0038] This equation shows that the contribution of specific ^^^^, ^^^^, and ^^^^ can be represented by the ratio of the partial sum of the given ^^^^, ^^^^, and ^^^^ to the total hazard. Note that ^^^^ and ^^^^ are constant. Thus, the hazard deaggregation, ^^^^(^^^^, ^^^^, ^^^^|^^^^ > ^^^^), is proportional to ^^^^(^^^^ > ^^^^|^^^^, ^^^^, ^^^^)^^^^(^^^^, ^^^^, ^^^^): ^^^^(^^^^, ^^^^, ^^^^|^^^^ > ^^^^) ∝ ^^^^(^^^^ > ^^^^|^^^^, ^^^^, ^^^^)^^^^(^^^^, ^^^^, ^^^^)

[0039] The Riemann summation This method computes PSHA curves by summing the areas of partitioned (^^^^, ^^^^, ^^^^) cuboids. The Riemann summation for Eq. (1) can be expressed as:^^^^^^^^^^^^^^^^^^^^^^^^(10) (^^^^ > ^^^^) = ^^^^^^^^(^^^^ > ^^^^|^^^^^^^^ , ^^^^^^^^ , ^^^^^^^^)^^^^^^^^,^^^^,ℰ (^^^^^^^^ , ^^^^^^^^ , ^^^^^^^^)Δ^^^^Δ^^^^Δ^^^^^^^^=1 ^^^^=1 ^^^^=1 TABLE 1: Comparison of Time Complexity of Various PSHA Algorithmswhere Δ^^^^, Δ^^^^, and Δ^^^^ are grid step size for Riemann sum, ^^^^^^^^, ^^^^^^^^, and ^^^^^^^^are the total number of grids, satisfying ^^^^^^^^Δ^^^^ = ^^^^max- ^^^^^^^^^^^^^^^^, where ^^^^ = {^^^^, ^^^^, ^^^^}. Note that^^^^^^^^, ^^^^^^^^)Δ^^^^Δ^^^^Δ^^^^ is equivalent to the probability at (^^^^, ^^^^, ^^^^) such that |^^^^ − ^^^^^^^^| < Δ^^^^∕2, |^^^^ − ^^^^^^^^| < Δ^^^^∕2, and |^^^^ − ^^^^^^^^| < Δ^^^^∕2. The accuracy of the Riemann summation depends on the grid size. Utilizing finer grids enhances the accuracy of the summation. However, the computation time is inversely proportional to the grid step size and thus proportional to the number of grids. For three-dimensional PSHA summation, the computation time scales with ^^^^^^^^ × ^^^^^^^^ × ^^^^^^^^ (Table 1). Notably, since the distance PDF (^^^^^^^^(^^^^)) cannot be analytically determined in practice, integration often extends over latitude (^^^^), longitude (^^^^), and depth (^^^^), increasing the dimensions from three (^^^^, ^^^^, ^^^^) to five (^^^^, ^^^^, ^^^^, ^^^^, ^^^^). Therefore, Riemann summation for seismic hazard becomes even more computationally intensive, a phenomenon known as "the curse of dimensionality" due to the exponential increase in computation time with the number of dimensions.

[0040] A Monte Carlo (MC) model generally uses random sampling of data to model the probability of different outcomes in a process that involves uncertainty. In PSHA, the MC method simulates many synthetic earthquakes and calculates PSHA by assessing the frequency with which ground motion intensities exceed a certain threshold. MC PSHA is computed aswhere ^^^^^^^^ denotes the simulated ground motions, and ^^^^ is the total number of samples. ^^^^^^^^, ^^^^^^^^, ℰ^^^^ are random samples from^^^^, ^^^^), ^^^^ represents the equivalent catalog duration equal to ^^^^∕^^^^. MC PSHA is unbiased because the expectation of ^̂^^^(X>a) is the same as ^^^^:=^^^^(^^^^ > ^^^^)(12) The variance of ^̂^^^, the variability of each MC estimates, can be obtained as follows: ^^^^^^^^ − ^^^^2(13) ^^^^^^^^^^^^�^̂^^^� =^^^^ Note that VAR[̂^^^^] is always positive since ^^^^ ≥ ^^^^. Then, the standard deviation of the estimate, ̂^^^^, can be obtained as: ^^^^^^^^ − ^^^^2(14)�^^^^ which provides the absolute uncertainty about MC PSHA estimates. However, one does not want to fully rely on ^^^^[^̂^^^] to compare different exceedance probabilities. For example, suppose someone is interested in two different exceedance frequencies, ^^^^1= 10−1 / yr and ^^^^2= 10−2 / yr when ^^^^ = 1 / yr and ^^^^ = 100. Then, the variance (uncertainty) of the two MC estimates, ^̂^^^1 and ^̂^^^2, are ^^^^1 = 3 × 10−2and ^^^^2 ∼ 1 × 10−2, respectively (Eq. (14). Here, someone might argue that uncertainty of ^̂^^^1 MC estimate is larger than that of ^̂^^^2 because ^^^^1 is greater than ^^^^2. However, the uncertainty of ^^^^1 should be considered smaller than that of ^^^^2 considering the target true value of each estimate. That is, 0.1 ± 0.03 (̂^^^^1) is better estimate than 0.01 ±0.01 (^̂^^^2). Instead, it is better to use the coefficient of variation, COV, the standard deviation normalized by its mean ^^^^[^̂^^^] ^^^^^^^^^^^^ = =�^^^^[^̂^^^] ^^^^^^^^ (15)

[0041] Smaller COV indicates that the PSHA estimate is more accurate relative to its actual value, e.g., 5% COV means that the hazard estimates are within ±5% of its true value with a 68% probability. COV is inversely proportional to the square root of the number of MC samples (^^^^), indicating that increasing ^^^^ naturally improves the estimate’s accuracy.

[0042] The target exceedance frequency, ^^^^, also affects the MC accuracy. The lower ^^^^ leads to poor accuracy with fixed ^^^^ and ^^^^. It is intuitively reasonable that sufficiently long earthquake catalog is required to accurately estimate the event from long return period. It is also notable that the variance of MC estimate is proportional to the square root of earthquake occurrence rate, ^^^^, indicating that MC PSHA is more challenging task in the region with higher earthquake activity. This can also be explained intuitively. Regions with high seismic activity experience more earthquakes per year than less active areas, therefore, when calculating the ground motion exceedance probability over the same period, more earthquakes need to be considered in these active regions.

[0043] An advantage of MC PSHA is that the computational time is determined by the number of samples, ^^^^, circumventing the dimensionality issue inherent in Riemann summation (Table 1). In other words, the computation time of MC simulation is not influenced by the number of grids (^^^^^^^^, ^^^^^^^^, and ^^^^^^^^ in Eq. (10)). Therefore, one could utilize fine joint probability mass functions for more precise hazard estimation without an increase in computational burden. If closed-form probability density functions are available, PSHA can be implemented without approximation by discretization. Also, the computation time of MC is independent of the number of ground motions of interest, ^^^^^^^^. The embodiments can utilize the generated synthetic ground motion catalog to tally the exceedance events for all the ground motions of interest, though one still should repeat the computation along with the number of MC samples.

[0044] Conventional MC faces extreme computational challenges for low probabilities because the number of samples required to achieve low COVs dramatically increases. One can compute the required number of samples by rearranging Eq. (15):1^^^^ − ^^^^(16) ^^^^ ~2× (^^^^^^^^^^^^) ^^^^ If ^^^^ >> ^^^^, the Eq. (16) can be approximated as:To demonstrate conventional MC’s extreme computational demands, the embodiments calculate ^^^^ for a target annual exceedance frequency, ^^^^, of 10−4per year and an annual earthquake rate ^^^^ = 1 / yr, which is typical values in PSHA practice. From Eq. (17), the process needs ^^^^ ∼ 108to achieve COV = 1%.

[0045] Importance Sampling (IS) is a generalization of MC theory that estimates an expectation by approximating it with a weighted average of random samples from another distribution. The embodiments use the IS distribution to sample rare events with a higher likelihood than conventional MC and then correct their frequency through weights, significantly reducing the number of samples to compute low probabilities. Consider a random variable ^^^^ that follows a probability function, ^^^^^^^^(^^^^). The expected value of a function ^^^^(^^^^), denoted as ^^^^, is defined byThe MC estimate is:where ^^^^^^^^is sampled from ^^^^^^^^(^^^^).

[0046] By introducing an arbitrary probability function, ^^^^^^^^(^^^^), Eq. (18) can be equivalently expressed as:The process restricts the integration range in Eq. (20) where ^^^^^^^^(^^^^) ≠ 0 because ^^^^ such that ^^^^^^^^(^^^^) = 0 does not contribute to the integration. Then, ^^^^^^^^(^^^^) can be any distribution with nonzero density in the integration range. Eq. (20) provides important implications in MC estimation. The embodiments obtain the solution to Eq. (18) by estimating the expected value of ^^^^(^^^^)^^^^^^^^(^^^^)∕^^^^^^^^(^^^^) where ^^^^ follows the distribution ^^^^^^^^(^^^^). This approach is highly useful innumerical integration, especially when sampling from ^^^^^^^^(^^^^) is challenging or the population of ^^^^^^^^(^^^^) is extremely low in the region of importance, e.g., when ^^^^(^^^^) is large.

[0047] The IS MC estimate is:where ^^^^^^^^s are sampled from ^^^^^^^^(^^^^), the proposed (or new) IS sampling density function. If ^^^^^^^^(^^^^) is equal to the original distribution, ^^^^^^^^(^^^^), Eq. (21) simplifies to the conventional Monte-Carlo (Eq. (19)). The ratio ^^^^^^^^(^^^^)∕^^^^^^^^(^^^^), known as the importance weight (^^^^^^^^), adjusts for the change in sampling distribution.

[0048] Researchers have applied IS to PSHA using different sampling functions. For IS PSHA, Eq. (1) can be reformulated by introducing a new sampling joint density function,The IS MC estimator of equation (22) isThe mean of ISAlso, the variance of the IS PSHA estimates with respect to the true ^^^^ is given byHence, the COV of IS PSHA estimate can be expressed as^^^^ (^^^^^(^^^^^^^^,^^^^^^^^,ℰ^^^)2 2^^^^^^^^^^^>^^^^|^^^^^^^^^^^^^^^^ℰ^^^^)^^^^,^^^^,ℰ ^^^^^ (^^^^^^^^,^^^^^,ℰ )�−^^^^2�^^^^^^^^^^ =�^^^^ ^^^^,^^^^,ℰ ^^^ ^^^^^^(25) ^^^^^^^^2From Eq. (24), the embodiments specifically choose a new sampling density ^^^^∗^^^^,^^^^,ℰthat makes VAR[̂^̂^^^]=0,Using ^^^^∗^^^^,^^^^,ℰ, one could compute the true hazard, ^^^^, with only one MC sample because IS MC is unbiased (Eq. (23) and q*the variance is zero. As used here ^^^^∗means the optimal IS density for PSHA calculation.

[0049] One should note that Eq. (26) is exactly the same as Eq. (9), indicating that the optimal sampling density, ^^^^∗, is identical to the hazard deaggregation. That is, if one finds ^^^^∗, the embodiments are able to not only dramatically enhance the computational efficiency of seismic hazard estimation but also obtain the hazard deaggregation distributions as a by- product.

[0050] In fact, one can also see this profound relationship between the hazard and deaggregation estimates by rearranging Eq. (9):

[0051] Because this identity holds for any values of (^^^^, ^^^^, ^^^^), one can obtain the hazard at ground motion level ^^^^ with any (^^^^, ^^^^, ^^^^) triplet if we know ^^^^(^^^^, ^^^^, ^^^^|^^^^ > ^^^^). The term ^^^^, ^^^^)∕^^^^(^^^^, ^^^^, ^^^^|^^^^ > ^^^^) can be interpreted as the importance weight of IS, and ^^^^(^^^^, ^^^^|^^^^ > ^^^^), the hazard deaggregation, is the optimal sampling density.

[0052] The deaggregation can be more computationally expensive than the hazard. In Riemann sum, one must save each sum element in memory and allocate those elements into appropriate deaggregation bins. Also, in conventional MC, one should save the long synthetic ground motion catalog with the corresponding (^^^^, ^^^^, ^^^^) triplet to allocate those into the proper bins. These operations necessitate significant computational memory and time. Therefore, the property that the optimal sampling density resembles the hazard deaggregation can be considered a big benefit for hazard analysts. IS MC with different densities other than ^^^^∗could also improve the computational efficiency compared to the conventional MC; however, it cannot give any information on the hazard deaggregation.

[0053] Though the use of optimal sampling density is a large benefit in computation of seismic hazard and hazard deaggregation, obtaining it is not a trivial because ^^^^(^^^^^^^^ > ^^^^|^^^^^^^^ , ^^^^^^^^ , ^^^^^^^^) and ^^^^ in Eq. (26) are unknowns before calculation of PSHA. Thus, the embodiments propose a new PSHA computation method to find it.

[0054] In “adaptive” importance sampling (AIS), the embodiments iteratively train the IS density to find the optimal one by exploring important regions to compute I(Xi> ^^^^ | Mi, Ri, Ei) and λ with a reduced number of MC samples. AIS must balance different factors in determining how many samples should be used to train optimal IS density. Fewer MC samples can make AIS fail as they will not allow effective exploration of these important regions effectively. On the other hand, many samples can impose large computational demands, even larger than those from conventional MC. In AIS, one should also consider computational costs are proportional to the number of iterations for convergence. Thus, it is important to select appropriate algorithms that converge fast.

[0055] VEGAS is a non-parametric adaptive importance sampling (AIS) algorithm that iteratively identifies the optimal proposal density, ^^^^∗. (Lepage, G. (1978). A new algorithm for adaptive multidimensional integration. Journal of Computational Physics 27(2), 192– 203.) The algorithm has been developed and widely used in computational physics, but it is currently applied to chemistry, astrophysics, finance, and medical statistics. The VEGAS algorithm refines the sampling strategy over several iterations, adjusting the sampling to use more samples in regions where the probabilities are higher.

[0056] The VEGAS algorithm is conceptually straightforward and is recognized for its rapid convergence, especially when the random variables involved are independent (Lepage, 1978, 2021). In PSHA, the variable ^^^^ is always considered an independent variable (Eq. (3). Additionally, the variables ^^^^ and ^^^^ are also treated as independent under point source assumption (Eq. (4). Notice, however, that when the finite-fault rupture model, in which the rupture dimension changes with magnitude, is adopted, ^^^^ and ^^^^ can be correlated, and the distribution of ^^^^ is conditional on the magnitude ^^^^.

[0057] In three-dimensional integration, which is the case of PSHA, VEGAS employs ^^^^3rectangular cuboids that are independently partitioned. The probability assigned to each cuboid and the total number of partitioned cuboids (^^^^3) are preserved across the iteration steps. However, the IS sampling density changes because the algorithm updates the cuboid’s size depending on its contribution to the integration. If a cuboid’s contribution is low, its size grows in the next step, lowering its probability density. Conversely, when a cuboid’scontribution is high, it shrinks in the following step, elevating the probability density. Ideally, when every cuboid’s contribution to the integration becomes identical, the embodiments find the proposed optimal sampling density, ^^^^∗, and the algorithm is terminated. The framework to find optimal IS densities for PSHA using VEGAS algorithm is explained with a simple point seismic source example in the following paragraphs.

[0058] First, the embodiments adopt an independently distributed joint probability function as IS density. Note that one could also include correlations in the IS density, but such approach would increase computational memory and time demands, e.g., the computational complexity increases exponentially with each added dimension, i.e., ^^^^(^^^^^^^^). In contrast, when one assumes independence, the computational complexity grows linearly with the number of dimensions, i.e., ^^^^(^^^^^^^^), making multi-dimensional integration in PSHA exceptionally efficient. Thus, ^^^^(^^^^, ^^^^, ^^^^) = ^^^^^^^^(^^^^)^^^^^^^^(^^^^)^^^^ℰ(^^^^) (27)Then, integration ranges for each variable – m, r, and ^^^^ – are divided into N grides with the same volume. This division is designed to generate cuboids of constant probability:

[0059] The number of grids, ^^^^, is chosen to be 50 as suggested by Lepage (1978) and because embodiments consider this value makes grids sufficiently fine to capture the actual distribution of the optimal sampling density ^^^^∗. The probability for each cuboid from Eq. (28) is set to 1∕^^^^3to ensure that the initial IS density function is uniformly distributed across the entire domain. As a result, the initial probability density ^^^^(0) for a specific cuboid isFIG.1 illustrates an example of initial partitioning when ^^^^=10. The gray dots in the FIG.1 shows an example of MC samples from this initial ^^^^0. Note that for easier visualization and understanding, FIG.1 illustrates an example in ^^^^-^^^^, the two-dimensional space, and not the three-dimensional ^^^^, ^^^^, and ^^^^ space.

[0060] The embodiments update ^^^^ using MC samples to ultimately make it converge to ^^^^∗(Eq. (26)).As mentioned earlier, the size of each cuboid is subject to change while the probability of each cuboid remains constant, i.e., probability density changes. This adjustment is facilitated through a "subdivision-and-restoration" process (Lepage, 1978). In this process, ^^^^th grid is subdivided into ^^^^^^^^sub-grids, with ^^^^^^^^ being proportional to the ^^^^th grid’s contribution to the overall integration, and restored to the original number, ^^^^, by merging ^^^^subgrid∕^^^^ consecutive subgrids, where ^^^^subgridis the total number of sub-grids. Thus, the number of subdivisions at ^^^^th grid is

[0061] ^^^^subgrid should be sufficiently larger than ^^^^ to iterate the IS density effectively, especially when the grids’ contributions to the hazard from the previous stage are highly heterogeneous. The embodiments chose ^^^^subgrid to be 10,000, which is 200 times greater than ^^^^ (=50). The second term of Eq. (29)-(31) represents the portion of each grid’s contribution to the hazard.

[0062] According to Lepage (1978), ^�^^^ in Eq. (29) is:where ^^^^(^^^^, ^^^^, ℰ; ^^^^) = ^^^^(^^^^ > ^^^^|^^^^, ^^^^, ℰ)^^^^^^^^,^^^^,ℰ(^^^^, ^^^^, ^^^^)

[0063] Note that ^^^^(^^^^, ^^^^, ℰ; ^^^^) is the integrand of PSHA. Intuitively, ^�^^^^^^^ can be considered asthe marginalization of the overall contribution of the magnitude dimension within the ^^^^th grid. Note that the summation of ^^^^2s in Eq. (32) is done over all the samples which fall within ^^^^^^^^−1 and ^^^^^^^^ and the division by ^^^^^^^^ ^^^^^^^^ can be considered an adjustment for unevenly distributed ^^^^ and ^^^^ samples to purely capture the magnitude contribution.

[0064] Similarly, ^�^^^^^^^ and ^�^^^^^^^ in Eq. (30) and Eq. (31) can be calculated as:^�^^^^^^^ = ∑ ^^^^2(^^^^,^^^^,ℰ;^^^^)^^^^^^^^−1<^^^^<^^^^^^^^^^^^^^^^(^^^^)^^^^ℰ(ℰ)^^^^2(^^^^,^^^^,ℰ;^^^^) ^�^^^^^^^ = �^^^^ <ℰ<^^^^^^^^^^(^^^^)^^^^^^^^(^^^^) ^^^^−1 ^^^^^^

[0065] In the early iterations, researchers have noticed numerical instabilities due to the relatively poor information on the integrand (Lepage, 1978, 2021). However, researchers have also found effective ways to mitigate it through smoothing (Eq. (33)) and damping (Eq. (34)).where ^^^^i is ^�^^^^^^^Δ^^^^^^^^∕ Σ^^^^ ^�^^^^^^^Δ^^^^^^^^ , and ^^^^ are ^^^^, ^^^^, or ^^^^ in equations 29-31.where ^^^^ is learning rate. The embodiments used ^^^^ as 1.0 as suggested by Lepage (2021). Through numerical experimenting, the inventors observed this value allowed most of the PSHA integration to converge within three iterations (as described later) without causing any numerical instabilities.

[0066] Finally, based on calculated ^^^^^^^^, ^^^^^^^^ and ^^^^^^^^ (Eq. (29)-(31)), the number of grids is restored to the original size, ^^^^, by merging ^^^^subgrid∕^^^^ consecutive subgrids. The restored grid is the updated density, ^^^^(1). FIG.1 illustrates this iteration. The sizes of the grids are shrunk at 5 <^^^^ < 6 and ^^^^ ∼ 2, indicating that the high contributions to the hazard on this range (Eq. (29)-Eq. (31)) in contrast to other less important regions, e.g., at ^^^^ > 6.5 and ^^^^ < 1. Note that the probability of each grid is preserved so that the number of MC samples inside each grid is almost the same regardless of the grid size.

[0067] The iterative process explained above is repeated until no further improvement is observed in the variance of the hazard estimator (Eq. (24). An example of the grid structure of the final proposed IS density is shown in FIG.1. The final proposed sampling density exhibits features which is expected based on intuition. First, the low contribution of ^^^^ less than 0 is understandable. This is because the simulated ground motion intensity, using Sadigh et al. (1997), never exceeds the target ground motion of 0.5 ^^^^ at any magnitude, even at ^^^^ =8.0, when ^^^^ is less than 0.1. In addition, the strong contribution of ^^^^ ∼ 5 and ^^^^ ∼ 2.2 to the integration makes sense because the simulated ground motion intensity exceeds the target ground motion of 0.5 g when ^^^^ and ^^^^ reach 5 and 2.2, respectively. The decreasing trend of contribution beyond ^^^^ ∼ 5 and ^^^^ ∼ 2.2 can be interpreted as the exponential decay in the probability of ^^^^ and ^^^^. The intuition from a simple PSHA example implies that the algorithm may work well for actual PSHA integration. TABLE 2: Seismic Source Information for Comprehensive Numerical Examples to Test the AIS PSHA Framework

[0068] The inventors tested the AIS PSHA embodiments on the comprehensive benchmark problem sets 1.11 and 2.1 from the PEER PSHA code verification project (Hale et al., 2018). For a thorough comparison, the inventors performed numerical computations using four different algorithms: Riemann sum, conventional Monte Carlo (MC), importance sampling (IS), and adaptive importance sampling (AIS) using the VEGAS algorithm. These computations were conducted across various seismic source settings: 1) an areal source (Area1), 2) a fault source (FaultA), and 3) a combination of an areal and two fault sources (Area1, FaultA, FaultB). These examples’ geometry and seismic activity parameters are detailed in FIG.2 and Table 2. FIG.3 presents the benchmark hazard curves for the examples. The ground motion model from Sadigh et al. (1997) was employed.

[0069] For Areal source 1, a circular-shaped areal source with a 100 km radius is considered (FIG.2). The earthquake activity rate of the source, ^^^^(^^^^ >^^^^min), is 0.0395 / year with ^^^^min and ^^^^max of 5.0 and 6.5, respectively. The earthquake occurrence model is doubly-truncatedexponential with Gutenberg-Richter ^^^^-value of 0.9. Seismogenic depth is 5 to 10 km from the surface, and a point source is assumed. Faults A and B are vertical fault sources (strike = 90◦E, dip = 90◦) with lengths of 50 and 85 km, respectively. They extend from the surface to 12 km depth. The earthquake activity models are chosen to be characteristic (Youngs and Coppersmith 1985) with the ^^^^-value of 0.9, and ^^^^minis set to be 5.0 for both fault sources. For Fault A, slip rate, ^^^^max, ^^^^^^^^ℎ^^^^^^^^are 1mm / year, 6.75, and 6.5; and for fault B, they are 2 mm / year, 7.0, and 6.75. In fault sources, ruptures are assumed to be floating inside the fault with the rupture dimensions following log10(^^^^) = ^^^^ – 4 (35) log10(^^^^) = 0.5^^^^ − 2.15(36)log10(^^^^) = 0.5^^^^ − 2.15 (37) where ^^^^ is rupture area in km2, ^^^^ is rupture width in km, ^^^^ is rupture length in km, and ^^^^ is earthquake magnitude. The probability of exceedance at eighteen ground motions (PGA), 0.001, 0.01, 0.05, 0.1, 0.15, 0.2, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 0.7, 0.8, 0.9, and 1.0 g, are estimated.

[0070] For Riemann summation, the embodiments followed the calculation procedure specified in the PEER PSHA code verification project Hale, C., N. Abrahamson, and Y. Bozorgnia (2018). Probabilistic Seismic Hazard Analysis Code Verification, PEER report 2018 / 03, Pacific Earthquake Engineering Research Center, Berkeley, CA.

[0071] Thus, the embodiments used magnitude and source-to-site distance step sizes of 0.01 and 0.1 km, respectively. Note that the number of grids of the source-to-site distance distribution (∼ 950) is significantly less than that of source location probability, the uniformly distributed probability of event location across the seismic source region, with the spatial and depth grid spacing of 0.5 km and 1 km ( ^^^^×1002 / 0.5×0.5) × 6 ≈ 750,000) suggested by one approach. This makes the Riemann sum computation ∼ 800 times more efficient. For the ground motion random variable (^^^^), the grid step size was set to 0.01, with minimum and maximum values of -6 and 6, respectively, beyond which the PSHA curve shows negligible change even at an annual exceedance probability of 10−8.

[0072] In the conventional MC approach, the embodiments generate a ground motion catalog based on the probability distributions of ^^^^, ^^^^, and ^^^^. From this catalog, the embodiments calculate the annual exceedance frequency by counting the instances where the groundmotion exceeds the specified threshold, as outlined in Eq. (11). This frequency is then converted to an annual exceedance probability using Eq. (6).

[0073] For IS, the embodiments adopt uniform IS densities across the entire integration range for magnitude, distance, and ^^^^. This choice is made due to the absence of prior information regarding which integration range significantly influence the hazard calculation.

[0074] For the embodiment of AIS PSHA with the VEGAS algorithm, the embodiments set the initial number of grids for ^^^^, ^^^^, and ^^^^ to 50 (Eq. (28)). The total number of subgrids, ^^^^subgrid, is chosen to be 10,000. The learning rate, ^^^^, is fixed at 1.0. Algorithm 1 presents the pseudocode. From line 1 to 5, the generation index ^^^^ is set to be zero, and the sampling density ^^^^ is set to be uniform. The main algorithm loop is from line 6 to 24. The loop is continued until there is no improvement in the coefficient of variation or there is no previous generation (line 6). At line 7, ^^^^^^^^, ^^^^^^^^, ^^^^^^^^ are sampled from the distribution ^^^^(^^^^). Then, the probability ^̂^^^ is estimated (line 8-10), and its coefficient of variation is also calculated (line 11). From line 12 to 23 the sampling function ^^^^(^^^^)is updated. The steps are repeated over ^^^^, ^^^^, and ^^^^ (line 12-15). From line 16 to 20, the contribution of each grid is calculated, at line 21, it is smoothed and dampened, and ^^^^(^^^^)is updated by subdivision depending on ^^^^^^^^and restore the number of grids to the original number, ^^^^ (line 22). Updated ^^^^(^^^^+1)is obtained by multiplying,^^^^r(^^^^+1), and ^^^^^^^^(^^^^+1)in line 23, and the while loop is ended by increase the generation index ^^^^. When main algorithm loop is terminated, it returns the hazard estimate ̂^^^^ and proposed optimal sampling density ^^^^ in line 25.

[0075] The inventors assessed the accuracy of conventional MC, IS, and AIS probability estimates through the standard deviation of the relative error with respect to the benchmark curve, calculated using the following formula:where ^^^^^^^^ is the relative error of the ^^^^th MC estimate, and ^^^^ is the total number of MC exceedance probability estimates. Since MC, IS MC, and VEGAS AIS estimates are unbiased (Eq. (12), (23)), the sum of ^^^^^^^^ when ^^^^ →∞ is theoretically zero. The relative error, ^^^^^^^^, is defined as:Algorithm 1 VEGAS adaptive importance sampling PSHA pseudocode [Parameters and Functions] ^^^^ Ground motion intensity of interest (e.g., ^^^^ = 0.1 g) ^^^^ Number of grids (e.g., ^^^^ =50) ^^^^subgrid Number of sub-grids (e.g., ^^^^subgrid =10,000) ^^^^^^^^ Number of MC samples (e.g., ^^^^^^^^ = 1,000) ⃗^^^^^^^^ Magnitude sample vector (⃗^^^^^^^^,^^^^ = ^^^^th element of ⃗X^^^^) ⃗^^^^^^^^ Distance sample vector ⃗^^^^^^^^ Ground motion random variable sample vector⃗^^^^ [ ⃗^^^^^^^^;⃗^^^^^^^^; ⃗^^^^^^^^]^^^^(⋅) Ground Motion Model ^^^^(⋅) Indicator function ^^^^^^^^(⋅) Original sampling distribution ^^^^(^^^^)(⋅) Proposed sampling distribution at ^^^^th iteration step ̂^^^^(^^^^)Hazard estimate at ^^^^th iteration step COV(^^^^)COV of the hazard estimate at ^^^^th iteration step ^^^^ pre-defined iteration stopping criteria (e.g., COV of 0.02) ^^^^ Index for partitioned grids (^^^^ = 1, 2,⋯, ^^^^) ^^^^ Index for samples (^^^^ = 1, 2,⋯, ^^^^^^^^) [Algorithm] 1: ^^^^ = 09: ⃗^^^^(^^^^) = ⃗^^^^(^^^^)∕^^^^(^^^^)( ⃗^^^^)10: ̂^^^^(^^^^)= Σ1^^^^^^^^^^^^1(^^^^)^^^^ ∕^^^^^^^^ 11: COV(^^^^) ←Eq. (25) 12: for ^^^^ in [^^^^, ^^^^, ^^^^] :13: if ^^^^ =^^^^ : (^^^^, ^^^^)←(^^^^, ^^^^) 14: if ^^^^ = ^^^^ : (^^^^, ^^^^)←(^^^^,^^^^) 15:16: for ^^^^ in {1, 2, ..., ^^^^} : 17: for ^^^^ in {1, 2, ..., ^^^^^^^^} : 18: 19:20: ^^^^^^^^ ←( ^^^^^^^^Δ^^^^^^^^ / Σ^^^^ ^^^^^^^^Δ^^^^^^^^) 21: ^^^^^^^^ ←smoothed, dampened ^^^^^^^^ (Eq. (33), 34) 22: ^^^^u(^^^^+1)←subdivision and restoration (Eq. (29)-31) 23: ^^^^(^^^^+1)←^^^^m(^^^^+1)^^^^r(^^^^+1)^^^^^^^^(^^^^+1)24: ^^^^←^^^^ + 1 25: returnwhere ^^^^^^^^is the exceedance probability estimated at ^^^^th MC estimate, and ^^^^^^^^^^^^^^^^is the benchmark probability. The analyses were conducted in a Python 3.11 on an Intel Core i7- 137002100 MHz processor with 64GB RAM.

[0076] The inventors considered the circular area source with the site at the circle’s center in FIG.2. For this example, the inventors show different computational performances to conduct PSHA in FIG.4. In FIG.4, conventional MC is shown in squares, IS MC in triangles, and AIS MC in circles. The numerical experiments show that computational times and standard deviations (“accuracy”) are linearly correlated in logarithmic scale, in agreement with the theory, because the required number of MC samples (linearly proportional to the computational time) is inversely proportional to the square of the standard deviation (Eq. (16)).

[0077] While the computation time and standard deviation vary depending on the target ground motion of interest, the results show the AIS PSHA generally outperforms the other numerical techniques. AIS PSHA becomes increasingly efficient for larger ground motions. At a low target ground motion of 0.05 g, the computation time to achieve a 2 % standard deviation is 0.01 seconds for AIS, while it takes 0.05 and 0.09 seconds for conventional MC and IS estimates, respectively, indicating AIS is 4.8 and 8.5 times faster. However, computational efficiency becomes extreme in high-ground motions. At 1.0 g, the computation time to achieve 2 % standard deviation is estimated to be 0.02 seconds for AIS, while it takes127 and 1.4 seconds for conventional MC and IS estimates. In this case, AIS is 7,800 and 70 times faster than conventional MC and IS, respectively. Also, note that for 2 % standard deviation case, AIS is > 105faster than Riemann sum. The inventors also note that AIS outperforms IS in all ground motion ranges by a factor of 8 to 70 to achieve a 2 % standard deviation. However, this is not always guaranteed because AIS takes ^^^^ times more computation time than IS with the same ^^^^ due to the ^^^^ iterations to find the optimal IS density. The findings imply AIS PSHA with the VEGAS algorithm can find (close to) optimal IS density quickly.

[0078] It is also noteworthy that AIS’s computational efficiency is quite similar across different target ground motions, while conventional MC’s efficiency decreases sharply for higher ground motions shown in FIG.5. In fact, if the hazard curve is exponential (Marzocchi and Jordan, 2017), one can show conventional MC decreases its efficiency also exponentially (Eq. (15)). For example, the error in the hazard curve for 1.0g will increase to 100% if it initially was 1 % for 0.001 g with exceedance frequency of 0.9 / year (when ^^^^ = 1 / year). From the numerical experiments, the inventors found errors grow from 0 to 5.60 % with ^^^^convMC = 107MC samples for these ground motions. In contrast, in case of AIS with a similar computation time, the inventors found AIS had errors in the same range varying between 0.13 and 0.64 %, for this quite different ground motion levels FIG.5.

[0079] This finding is key for PSHA as the computational bottleneck is at the highest ground motion intensity. Using conventional MC PSHA, the hazard analyst has no choice but to largely increase the number of samples to estimate hazard accurately at high ground motions even though such a large number would not be necessary for low ground motions. This makes the conventional MC highly inefficient, and consequently, the efficiency of the conventional MC for lower ground motions cannot be considered a real advantage for PSHA. AIS PSHA overcomes this problem by adopting different optimal densities at different ground motion intensities, making the computational burden almost flat for any ground motion intensity, as shown in FIGs 6A-6D.

[0080] In terms of the accuracy of AIS PSHA, the PEER PSHA verification project suggests a strict acceptable error range of 5% for reliable PSHA computation codes. AIS PSHA achieves it with only N ∼ 10,000 per ground motion, as shown in FIGs.6A-6D.

[0081] Another key advantage of the embodiments of AIS PSHA is the co-production of deaggregation curves at no extra computational cost. From Eq. (26), optimal IS densities are theoretically equivalent to hazard deaggregation distributions. The results show AIS can findclose-to-optimal IS densities and thus closely match hazard deaggregation curves. The inventors found the benchmark marginal distributions of hazard deaggregation obtained from Riemann sum closely match the iterated IS density from AIS PSHA (FIG.7). The embodiments used the Kolmogorov-Smirnov (K-S) ^^^^ statistic (Kolmogorov, 1933) to quantify their similarities. K-S ^^^^ statistic measures the maximum difference between two cumulative distribution functions (CDF). If two CDF are identical, ^^^^ is zero, and its maximum possible value is one. ^^^^ close to zero indicates that two probability distributions are similar. The process calculated ^^^^ for ^^^^, ^^^^, and ^^^^ at all the ground motion intensities as shown in FIG.8. The process finds the maximum ^^^^ values for ^^^^, ^^^^, ^^^^ were 0.032, 0.113, and 0.092, respectively, and the minimum values were 0.019, 0.026, and 0.017, indicating a strong resemblance between the two distributions.

[0082] The method of the embodiments also estimated the differences of mean values from proposed ^^^^∗ and Riemann sum hazard deaggregation as shown in FIG.8. Note that the deaggregation distributions and optimal IS densities vary for different ground motion levels. Thus, their mean values also vary. The method of the embodiments found the maximum relative differences were 2.5 %, 22.6 %, and 18.5 % for ^^^^, ^^^^, ^^^^, and the minimum differences were 1.7 %, 1.0 %, and 4.3 %, respectively. Also, the maximum absolute differences were 0.14, 3.9 km, and 0.17 for ^^^^, ^^^^, ^^^^, and the minimum differences were 0.09, 0.7 km, and 0.0009. The maximum relative difference in distance (^^^^) appears for deaggregation distributions at 0.25 g, where the mean distance obtained from hazard deaggregation is 14.3 km and that from ^^^^∗ is 17.5 km. Given that users are typically interested in distance ranges on the order of tens of kilometers (e.g., 0-15 km, 15-25 km, 25-50 km, etc.) rather than a single value (U. S. Nuclear Regulatory Commission, 2007), this difference is not crucial in determining the controlling earthquake for critical infrastructures. In addition, IS densities still show small K-S ^^^^ statistics and their mode almost matches each other, as shown in FIG. 7.

[0083] As an experiment, the inventors considered a 50 km-length vertical fault 25 km away from the site as shown in FIG.2. The method adopted a finite-dimension rupture model, which results in a distance distribution dependent on magnitude. The results indicated this dependency can diminish the performance of VEGAS AIS because the VEGAS algorithm assumes the independently distributed random variables. For example, for a ground motion intensity of 0.05 g, AIS is slower than conventional MC, as shown in FIG.9, while AIS outperformed conventional MC at the same ground motion intensity when point sourceassumption was made, as seen earlier in FIG.4. In FIG.9, the conventional MC is shown in square, IS MC in triangles, and AIS MC in circles. However, as the ground motion intensity increases, the computational gap between the two methods becomes smaller rapidly and closes at 0.2g. For higher ground 413 motions, AIS outperforms conventional MC. At 1.0 g, for a 5% standard deviation, AIS, IS, and conventional MC take 0.08, 3.7, and 166 seconds, respectively, i.e., AIS MC is 48 and 2,162 times faster.

[0084] The embodiments also present the accuracy of AIS PSHA in comparison to the true hazard calculated by Riemann summation, shown in FIGs.10A-10D. Implementing the methods, the inventors observed AIS PSHA estimates approximate the true hazard curve within an acceptable range (5 % error) when ^^^^ is greater than ∼ 50,000. Also, note that though the embodiments assumed the independently distributed optimal sampling density, AIS PSHA still gives an unbiased hazard estimate that will converge to the true hazard with a sufficiently large number of samples due to the nature of the importance sampling, shown in FIG.10D.

[0085] The results also compared hazard deaggregation in the black line and the iterated IS density in the gray line in FIG.11 and showed they closely match each other even though the distance distribution depends on the magnitude in this case. It is also noteworthy the iterated IS density can even reproduce complex densities with discontinuities like the large jump within magnitude distribution (for ^^^^ =6.25) due to the use of a characteristic earthquake occurrence model. The K-S D statistic and mean difference of the two distributions are also presented in FIG.12. The maximum values of D in ^^^^, ^^^^, and ^^^^ are 0.30, 0.68, and 0.13, respectively, and the minimum values are 0.02, 0.31, and 0.04. The inventors found that the largest discrepancies occur in the magnitude distribution, but errors can be considered negligible as the mean magnitude difference is within 6 % error. The inventors also noted considerable discrepancies in the distance distribution shape (see K-S D statistics) as curves with concentrated probabilities in narrow ranges are harder to estimate for AIS (FIG.11). However, the mean distance is still within a 5 % error range. The ^^^^ is generally in good agreement across all ground motion intensities. The results indicate higher differences at lower ground motion because the mean ^^^^ is close to zero. For example, the calculated mean ^^^^ at ground motion intensity of 0.01 g is 0.022 and 0.007, which is not a large difference in practice.

[0086] The maximum relative differences in mean values were found to be 6.2 %, 5.3 %, and 129 % for ^^^^, ^^^^, ^^^^, and the minimum differences were 0.006 %, 3.9 %, and 2.2 %,respectively. Note that the 129 % of ^^^^ case is corresponding the case where the mean ^^^^ is close to zero. The relative difference appears to be slightly higher than previous areal source example, however, the absolute difference is still remain to be significantly small, the maximum absolute differences were 0.36, 1.3 km, and 0.20 for ^^^^, ^^^^, ^^^^, and the minimum differences were 0.0004, 1.0 km, and 0.07.

[0087] In PSHA, one often has multiple seismic sources. The embodiments here considered one areal and two fault sources around the site to represent this case. This application posits a different mathematical problem than the previous two examples because one must introduce an additional variable to formulate AIS.

[0088] First, the probability of earthquake occurrence at ^^^^th seismic source can be defined as:where ^^^^^^^^ is the number of seismic sources (^^^^^^^^ = 3 in this example), and ^^^^^^^^ is the annual earthquake occurrence rate of ^^^^th seismic source. Because the discrete random variables cannot be used in AIS, the embodiments define a continuous random variable and its corresponding probability density function as^^^^^^^^(^^^^) is a piece- wise constant function where the heights are proportional to the corresponding sources’ earthquake occurrence rates. The embodiments introduce ^^^^^^^^(^^^^) into the PSHA integration and obtainwhich is the integral version of equations 5 and 1. Solving Eq. (39) posits computational challenges than the single source problems because ^^^^^^^^(^^^^) have the large jumps at ^^^^ = 1, 2, … , ^^^^^^^^ − 1. In addition, the embodiments introduce additional dependencies in the multi-source case because ^^^^ and ^^^^ depend on ^^^^, further diminishing the VEGAS algorithm’s effectiveness. Thus, the process tests two AIS PSHA approaches for this case: 1) full AIS approach utilizing the Eq. (39) and 2) partial AIS approach, which is the simple summation of single-source AIS PSHA curves.

[0089] The inventors compared the computational performance of MC, IS, and full and partial AIS shown in FIG.13, with conventional MC in light green, IS MC in yellow, partial AIS in bright green, and full AIS in red. Like the previous examples, conventional MC is faster for low ground motions, e.g., 0.05 g, but not high ones. For example, to achieve the ^^^^ of 2 % at a ground motion of 0.8 g, conventional 454 MC, IS, and full and partial AIS take 73.8, 10.0, 0.77, and 0.13 seconds, respectively, i.e., partial AIS is the most efficient algorithm, ∼ 583, 79, and 6 times faster than conventional MC, uniform IS, and full AIS, respectively.

[0090] As stated earlier, full AIS is less efficient than partial AIS because of the additional dependencies and jumps introduced by ^^^^^^^^(^^^^). FIGs.14A-14D and 15A-15D present the accuracy of both AIS approaches with different sample sizes. The contrast between the graph in FIG.14C, and the graph in FIG.15C, shows that partial AIS has smaller error than full AIS even with fewer samples (150,000 < 500,000). The results do not show the comparison for the partial AIS on ground motions greater than 0.8 g because fault B reaches a numerical instability due to its low exceedance probability (< 10−12 / yr), shown in FIGs.15A-15D.

[0091] FIG.16 shows a flowchart of the overall method of performing PSHA. This method will generally involve a computing device to implement the steps of the method, using a computing device having one or more processors configured to execute code to perform the method of the embodiments.

[0092] The PSHA method of the embodiments calculates the ground motion intensity exceedance probability (the mean hazard) at the site of interest given pre-specified ground motion intensity (a), such as 0.1 g; g = gravitational acceleration ~ 9.81 m / s2), and the hazard deaggregation, which is the probability distribution showing the relative contribution of each basic random variables of the PSHA to the mean hazard.

[0093] At 10 the method starts with the input of the seismic source model and the ground motion model. The seismic source model specifies the probabilistic distribution of how large (magnitude, ^^^^) and where (location or distance from the site, ^^^^) the earthquakes will occur. The ground motion model provides the probabilistic distribution on how large ground motion intensity (spectral acceleration) will be observed at the site of interest given the magnitude (^^^^) and location (^^^^) of the earthquake. The model is normally expressed as log-normal distribution with the parameters μ and σ. To simulate the random ground motion, standard normal random variable ε incorporated. The seismic source models and ground motion models are represented as model parameters, ^^^^, which are also defined as a probabilitydistribution depending on the likelihood that the parameter values are true. This is the third input of PSHA. At 12 the method constructs the joint (or conditional) probability distribution,^^^^^^^^(^^^^, ^^^^, ^^^^, ^^^^).

[0094] Using the adaptive importance sampling algorithm (e.g., VEGAS) of theembodiments the method obtains the optimal sampling density,^^^^, ^^^^, ^^^^), at 14 basedon the characteristic that ^^^^∗^^^^ (^^^^, ^^^^, ^^^^,is proportional to ^^^^(^^^^ > ^^^^|^^^^, ^^^^, ^^^^, ^^^^)^^^^^^^^(^^^^, ^^^^, ^^^^, ^^^^), where ^^^^ is the simulated ground motion intensity, ^^^^(^^^^ > ^^^^|^^^^, ^^^^, ^^^^, ^^^^) is theindicator function that takes 1 if ^^^^ > ^^^^, otherwise, 0.

[0095] At 16, the method computes, the mean hazard. lambda (λ), using Monte-Carlo simulation, ^^^^ is the annual earthquake occurrence, based on numerically obtained optimal sampling density obtained in 16. Also, the method computes the hazard deaggregation bymarginalizing ^^^^∗^^^^ (^^^^, ^^^^, ^^^^, ^^^^) with respect to ^^^^, ^^^^ and ^^^^.

[0096] At 20, the process ends with the output is the mean PSHA estimate, ^^^^, and the optimaldistribution which represents the PSHA hazard deaggregation, ^^^^(^^^^, ^^^^, ^^^^).

[0097] The embodiments involve a novel computational method for PSHA. The embodiments formulated the new AIS PSHA framework implementing the VEGAS algorithm. The above discussion compared the proposed method with existing computational frameworks: a) Riemann sum, which exhibits exponentially increasing computational cost as the grid interval becomes finer and high sensitivity to the grid spacing strategy; b) conventional MC, which require a substantial number of repetitive ground motion simulations, particularly for hazards with low exceedance probabilities; and c) importance sampling with uniform density, which lacks prior information on important regions in PSHA integration.

[0098] The results indicate that AIS PSHA outperformed all existing PSHA computational frameworks. AIS PSHA can dramatically reduce computational times by up to factors of > 105compared to traditional Riemann summation. It is also faster by a factor of 103compared to the conventional MC while maintaining the same level of accuracy. Additionally, it was also up to 70 times faster than IS with uniform sampling function, demonstrating the effective convergence of VEGAS algorithm in PSHA computation.

[0099] AIS PSHA requires similar computation time for any ground motion level, making its application to larger ground motions with low probability substantially more efficient than with conventional MC, i.e., easier extension of the hazard curve to larger ground motions. For the extension, using conventional MC, one should inevitably generate larger number ofsamples which is not necessary in estimating the hazard at lower ground motion ranges. However, AIS PSHA framework trains the sampling density optimized to each ground motion level, minimizing the number of required samples regardless of the target ground motion level.

[0100] AIS PSHA finds approximated deaggregation curves at no extra computational cost based on theoretical insights showing that optimal IS densities are equivalent to deaggregation distributions. The above discussion shows empirically that the hazard deaggregation and iterated IS densities from AIS PSHA are fairly similar by comparing the statistical properties of the two distributions, e.g., K-S D statistics < 0.113 and mean values differences of <4.3 %.

[0101] The finite rupture model, which presents the dependency of the magnitude and distance distribution, can also be calculated using the VEGAS AIS algorithm while maintaining computational efficiency and the resemblance of the proposed density to hazard deaggregation.

[0102] In the case of combined seismic sources, the embodiments include two strategies for multiple source PSHA computation: 1) incorporating source probability into AIS framework and 2) simple summation of individually computed AIS PSHA curves. Both strategies outperformed the traditional methods up to by a factor of ∼ 580. The simple summation strategy is found to be more efficient than incorporating the source probability into the AIS framework by a factor of 6, implying a limitation of VEGAS algorithm in computation of integration involving dependent variables.

[0103] AIS can be applied to any PSHA computation. Specifically, for large-scale PSHA projects that involve numerous logic tree branches, the suggested method can be highly useful in reducing extensive computational demands.

[0104] BEGIN LOGIC TREE DISCUSSION

[0105] Another embodiment addresses issues with input models to PSHA having inherent uncertainties due to limitations in current scientific understanding. These model-related uncertainties are termed epistemic uncertainties. To account for them, parameter distributions are defined and multiple alternative hazard curves are generated. For instance, hazard curves are evaluated across a range of seismic source parameters such as the maximum magnitude, Gutenberg-Richter b-value, earthquake rate, and across alternative ground motion models. The mean hazard curve is then computed by averaging these alternatives, while fractilehazard curves (e.g., 16th, 50th, and 84th percentiles) are used to represent the uncertainty range. These mean and fractile hazard curves are critical outputs in modern PSHA practice.

[0106] Numerical evaluation of the PSHA mean and fractiles is traditionally achieved via double-nested integration The inner integration evaluates aleatory uncertainty, the individual hazard curves, typically using Riemann sums or Monte Carlo (MC) simulations. The outer integration propagates epistemic uncertainty in model parameters such as the maximum magnitude, b-value, and ground motion model parameters. The latter is commonly implemented via a logic-tree framework, where branches represent alternative models or parameters and are weighted based on expert judgment.

[0107] Despite its widespread adoption, the logic-tree approach presents limitations in both computational efficiency and accuracy. From a computational standpoint, the nested integration process is highly demanding. Two primary factors contribute to the computational burden: 1) the time-intensive evaluation of individual hazard curves, especially at low exceedance probabilities, and 2) the exponential growth in the number of individual hazard evaluations required as the number of epistemic variables increases. This "curse of dimensionality" often forces analysts to simplify the logic tree by reducing the number of branches, which, in turn, necessitates additional sensitivity analyses, 36 further increasing computational costs.

[0108] In terms of accuracy, the logic-tree method approximates continuous epistemic uncertainty distributions using discrete representations. This discretization can introduce bias in the resulting hazard curves, especially when the number of discrete points per variable is limited, typically to three to five. Furthermore, the process is subjective, as different experts may pro-40 duce different discretizations for the same underlying distribution. Although strategies exist to reduce discretization bias, they cannot eliminate it entirely. Increasing the number of discretization points could improve accuracy, but this significantly increases computational cost. Specifically, the complexity scales exponentially as O(nd), where n is the number of discretization points per variable and d is the number of independent epistemic variables. For instance, increasing the number of points from three to five across ten variables results in a 165-fold increase in computational cost (510 / 310), underscoring the interdependence between accuracy and efficiency in PSHA.

[0109] To mitigate these challenges, several studies have explored Monte Carlo-based importance sampling (IS) methods. IS improves efficiency by focusing sampling efforts on regions that contribute most to the integration result, thereby reducing the number of requiredsimulations. However, identifying these important regions in PSHA is nontrivial, motivating the use of adaptive importance sampling algorithms that iteratively refine the proposal distribution. The embodiments above regarding sampling involve an adaptive IS framework that closely approximates the optimal sampling density and significantly accelerates individual hazard and disaggregation computations. However, this approach assumes independence among variables, an assumption invalidated in the presence of epistemic dependencies, rendering it unsuitable for 54 computing mean and fractile hazards.

[0110] One proposal the uses polynomial chaos expansion (PCE) to propagate epistemic uncertainty in ground motion models. This method effectively reduced computational costs, particularly for non-ergodic models, by providing surrogate models of ground motion medians. However, its applicability is limited, as extending PCE to other 58 epistemic variables—such as seismic source parameters and ground motion standard deviations— requires closed-form expansions, 59 which are not always attainable.

[0111] The embodiments herein present a generalized computational framework for computing the mean and fractile hazard curves along with sensitivity analysis in PSHA using adaptive importance sampling. The method builds upon the adaptive IS framework developed and disclosed above, extending its applicability beyond individual hazard and disaggregation computations. The process begins by mathematically formulating PSHA to establish the groundwork for our approach. The bias and inefficiency introduced by the logic-tree framework is then analyzed, particularly due to its discretization of continuous distributions. The embodiments then introduce a generalized adaptive IS method and incorporate it into the PSHA formulation to develop a novel framework based on Gaussian Population Monte Carlo Adaptive Importance Sampling (G-PMC AIS). The proposed approach efficiently estimates both mean and fractile hazards, while enabling sensitivity analysis. The discussion demonstrates its computational advantages through numerical examples that incorporate commonly used epistemic uncertainties, including the uncertainties in both seismic source and ground motion models.

[0112] Given the seismic source and ground motion models, the individual hazard for a single source can be expressed as:where θ denotes the model parameters associated with epistemic uncertainty, such as the Gutenberg–Richter b-value, mmax, and alternative ground motion models. The variables (m, r, ε) represent magnitude, distance, and the ground motion residual, respectively, encapsulatingthe aleatory uncertainty in hazard. The function^^^^|^^^^) denotes the conditionalprobability density function (PDF) of the aleatory variables given θ. The term ν(θ) represents the earthquake occurrence rate under θ, and u(m, r, ε; θ ) is the simulated ground motion intensity. The indicator function I(·) evaluates to 1 if u(·) exceeds the ground motion level a.

[0113] The aleatory (random) variables (m, r, ε) can be collectively denoted as φ, leading to the simplified expression:

[0114] The mean hazard is defined as the expectation of the individual hazard over the epistemic uncertainty variables θ:

[0115] Here, λ (a) represents the mean hazard, and ^^^Θ^ (^^^^) is the PDF of the epistemicuncertainty variables. Note that equation Eq. (3) corresponds to the discrete form used inlogic-tree-based hazard computation, λ (a) =^^^^^^^^(^^^^^^^^), where N^is the total number of logic-tree end branches, anddenotes the weight associated with the i-th branch.

[0116] The distribution of individual hazard with respect to θ, denoted λ(a|θ), represents the range of possible hazard values. Its cumulative distribution function (CDF) is given by: ^^^^^^^^(^^^^) = ^^^^[^^^^(^^^^|^^^^) ≤ ^^^^] (4)(43)In practical applications, the probability on the right-hand side of Eq. (4) is computed empirically by cumulatively summing the weights in ascending order of individual hazard values.

[0117] The p-th fractile of the hazard, denoted λp(a), is obtained from the CDF as the smallest value for which the CDF reaches or exceeds p%: ^^^^^^^^(^^^^) = inf�^^^^:^^^^^^^^(^^^^)(5)(44)where λp(a) denotes p-th fractile of the hazard, ^^^^^^^^(^^^^)is the CDF of the individual hazard, and inf{·} identifies the minimum hazard value satisfying the fractile condition. The accuracy ofλp(a) depends on the precision of the estimated ^^^^^^^^(^^^^).

[0118] The discussion now turns to examining the issues of accuracy and computational efficiency associated with the traditional approach to propagating epistemic uncertainty in PSHA: the logic tree. The embodiments address the case where the continuous distributions of the epistemic uncertainty variables are approximated using discrete distributions through logic trees.

[0119] The embodiments show the potential bias introduced by discretizing a continuous distribution through a simple numerical example, shown in FIG.17A and 17B. Consider a point source located 35 km from the site that generates a single characteristic earthquake. For illustration, the simple base ground motion model (GMM) by Sadigh et al. (1997) is used. Three epistemic uncertainty variables are incorporated: the characteristic magnitude, median ground motion, and standard deviation of the GMM. The underlying distribution for the characteristic magnitude is assumed to be normally distributed with a mean of 7.0 and a standard deviation of 0.2, ^^^^ (7.0, 0.2). The natural logarithm of the median and the standard deviation (ln) of ground motion are modeled as normally distributed deviations from the base GMM of Sadigh et al. (1997), both with standard deviations of 0.05, i.e., ^^^^(0, 0.05) and ^^^^ (0, 0.05), respectively. The distributions are not truncated to upper or lower limits since the ranges yielding unrealistic values, e.g., negative standard deviation, do not affect the result at the hazard range considered in this example. The annual occurrence rate of the characteristic earthquake is set to 0.01.

[0120] To construct the logic tree, the embodiments consider three discretization strategies (Table 3). The Extended Pearson-Tukey three-point approximation of Keefer and Bodily (1983) (denoted KB83), the three-point approximation used in the 2023 USGS seismic hazard model (Petersen et al., 2024) (denoted Pea24), and the five-point approximation proposed by Miller and Rice (1983) (denoted 111 MR83). These approximation schemes have been widely used in many practical applications, and they aim to match the first few statistical moments of the original distribution.

[0121] The benchmark solution is obtained using continuous-distribution-based Monte Carlo(MC) simulation with N = 5 × 108 samples, which ensures a coefficient of variation (COV) ≤1% down to a hazard level of 2 × 10−7 / yr.

[0122] The mean, 16th, and 84th fractile hazards from each approximation scheme, along with their relative errors compared to the benchmark, are shown in FIGs.18A-18C. As expected, increasing the number of discrete points improves accuracy. For the five-point MR83 approximation, errors remain within 5.28 % in ground motion intensity up to 0.5 g for the mean, 16th, and 84th fractile hazards. However, a tendency to diverge from the benchmark at higher ground motions is still observed.

[0123] The KB83 approximation yields mean hazard estimates comparable in accuracy to the five-point approximation, despite requiring only 22 % of the computational cost (= 33 / 53). It also produces reasonably accurate 84th fractile estimates, with errors under 15 %, and does not show significant divergence at higher ground motion levels within the range considered. However, for the 16th fractile, KB83 significantly underestimates hazard values beyond 0.1 g (exceedance rate of 6.5 × 10−3 / yr). Although this exceedance level corresponds to a relatively high probability event (probability of exceedance without accounting for earthquake rate is 0.65), the KB83 approximation still deviates substantially, indicating difficulties in accurately capturing the lower tail of the hazard distribution. Note that KB83 underestimated the 16th fractile, suggesting an exaggerated range in the hazard curve, driven by a low bias in the lower fractiles.

[0124] Among the three, Pea24 performed the worst. The mean hazard deviates from the benchmark by over 5% at ground motion levels exceeding 0.2 g (hazard level of 1.1 × 10−3 / yr). The 16th and 84th fractile estimates diverge even earlier, at ground motion levels of 0.1 g, corresponding to hazard levels of 1.8 × 10−3 / yr and 3.9 ×10−3 / yr, respectively. Notably, Pea24 overestimates the 16th fractile and underestimates the 84th fractile, suggesting that the scheme underrepresents the full variability of the hazard. Although both KB83 and Pea24 use the same number of branches, their differing weighting schemes lead to substantially different hazard estimates. At 0.5 g, the mean, 16th, and 84th fractile hazard estimates differ by factors of approximately 1.8, 2.7, and 1.8, respectively.

[0125] This example highlights that discretization in logic trees introduces bias into hazard estimates, especially at large ground motions (low exceedance probabilities) and in the lower fractile estimates. Any existing approximation scheme is not sufficient to accurately propagate the full distribution of the underlying continuous distribution, i.e., they do not fully capture the entire range of the hazard. Moreover, the example demonstrates that discretization strategies without statistical justification can result in significantly biased 140 hazard estimates even when they approximate the same underlying distribution.

[0126] In summary, this example reveals several limitations of logic-tree-based approximations in PSHA computation: First, the under lying continuous distributions of epistemic uncertainty variables can be significantly distorted when arbitrary weights are assigned, as shown in FIGs.17A-17B. Second, even the most statistically grounded discretizations cannot fully represent the underlying distributions, introducing errors in hazard estimates. This issue could be further exacerbated when a few epistemic variables dominate hazard sensitivity. In such cases, discretizing these key variables using only a few points causes an abrupt change in the hazard only when their values change, resulting in a step-like distribution of hazard that deviates significantly from the true fractile.

[0127] In the previous numerical example, the process only considered three epistemic uncertainty variables in the hazard calculations, requiring 27 and 125 hazard evaluations for the three- and five-point discretization schemes, respectively. However, practical PSHA applications involve more epistemic variables, resulting in substantial computational demands even after discretization. For instance, the total number of end branches in the 2023 USGS hazard maps for Hawaii is approximately 2.9 trillion, rendering direct computationinfeasible. Also, for less extreme cases, PSHA using UCERF3 results in ∼ 170,000 logic treeend branches, which is also required to be reduced to calculate the risk to a portfolio of buildings. These examples clearly illustrate the curse of dimensionality inherent in logic-tree- based propagation of uncertainty, where the computational cost 156 increases exponentially with the number of dimensions (O(nd )).

[0128] To demonstrate this effect, the embodiments present the estimated computation time for full hazard evaluations with increasing numbers of epistemic uncertainty variables, assuming an individual hazard computation time of five seconds. With ten epistemic variables and a five-point discretization, the total number of hazard evaluations becomes 510, making the computation time without parallelization approximately one year.

[0129] Another critical computational consideration—often overlooked in PSHA practice— is the number of marginal hazard evaluations required for accurate computation of each individual hazard curve. To ensure the reliability of the individual hazard estimate, a sufficient number of ground motion slices (in Riemann sum approaches) or samples (in Monte Carlo approaches) must be used, and the required number varies depending on the exceedance probability associated with the target ground motion level. Table 3: Approximate strategies commonly used in PSHA, testing in the numerical example for logic-tree

[0130] For example, consider a case incorporating uncertainty in the median ground motion model. Assume that three logic-tree branches are defined: a base case with median ground motion of 0.1 g, and upper and lower cases representing ±0.3 shifts in the logarithmic ground motion, corresponding to median ground motions of 0.075 g and 0.135 g, respectively. Assuming an aleatory variability of σ = 0.6 (natural logarithm), the exceedance probabilities at a target ground motion of 0.246 g are 15.9%, 6.7%, and 2.3% for the base, upper, and lower cases, respectively—a variation by a factor of 7. This discrepancy becomes even more pronounced at higher ground motion levels. For instance, at 0.45 g, the exceedance probabilities become 2.23%, 0.62%, 171 and 0.135%, resulting in a variation by a factor of 17. Such differences in exceedance probability directly influence the number of sufficient slices in Riemann sum or sufficient samples required in Monte Carlo simulation to achieve a consistent level of accuracy across branches. For example, the number of Monte-Carlo samples for the accurate calculation is inversely proportional to the correspondingexceedance probability (N ∼ 1 / λ ). Thus, to maintain equivalent accuracy, the lower-casebranch would require 17 times more samples than the upper-case.

[0131] One embodiment of Adaptive Importance Sampling was discussed above. Conventional Monte Carlo (MC) simulation comprises a general numerical method for estimating the expected value of a function. Consider a random variable X that follows a PDF fX (x). The expected value of a function g(x), denoted by S, is defined as: (6)(45)The conventional Monte Carlo estimate of S is given by.where Xi is the i-th random sample drawn from the distribution fX(x), and N is the total number of samples. As N increases, ^̂^^^^^^^^^^^converges to the true value S. The accuracy of the Monte Carlo estimate ^̂^^^^^^^^^^^can be quantified by the variance of the estimator, ^^^^^^^^^^^^�^̂^^^^^^^^^^^�:

[0132] The accuracy of the estimator is often represented by the standard deviation of the estimate divided by the mean, the coefficient of variation (COV):

[0133] As the COV decreases, the estimate becomes more robust. For example, a COV of 5% indicates that there is a 95% confidence (Z-score of 1.96) that the true value of S lies within ± 9.8% (= 1.96 × 5%) of the estimator ^̂^^^^^^^^^^^. In contrast, a COV of 1% implies a narrower range of ± 1.96% (= 1.96 × 1%) with the same level of confidence.

[0134] Importance Sampling (IS) is a generalization of the conventional Monte Carlo method. Eq. (6) can also be rewritten as:where qX(x) is the proposal density (or distribution), an arbitrary PDF that is non-zero over the domain of integration. The IS estimate of S is given by:where the random sample Xiis drawn from the proposal distribution qX(x), and N denotes the total number of samples. ^̂^^^^^^^^^^^is also an unbiased estimator, and the ratio fX(Xi) / qX(Xi) represents the importance weight, which adjusts for the change in the sampling density. The variance of the IS estimate is given by: (10)(49)The accuracy metric of SˆIS, COV, can also be expressed as:By properly selecting qX(x), COV [^̂^^^^^^^^^^^] can be reduced. Specifically, one can choose qX(x) to be:

[0135] This distribution ^^^^^∗^^^(^^^^) is referred to as the IS optimal sampling density since the variance and COV of the Importance Sampling estimate becomes zero (see Eq. (10)). This indicates that the S can, in theory, be estimated with a single random sample if one knows ^^^^^∗^^^(^^^^) at any x. However, identifying the optimal sampling density is nontrivial, as ^̂^^^^^^^^^^^is unknown and closed-form expressions for g(x) and fX(x) are not always available. To address this, the embodiments employ an algorithm that iteratively explores the region of importance, discussed above, Adaptive Importance Sampling (AIS), more particularly. Population Monte Carlo (PMC) AIS to apply it to PSHA.

[0136] Numerous AIS algorithms have been developed to iteratively explore regions of importance in numerical integration. Among them, Population Monte Carlo (PMC) is a straightforward AIS algorithm that utilizes importance weights to update the proposal distribution. Similar to other iterative algorithms, each iteration of the PMC AIS consists of three 212 fundamental steps aimed at approximating the optimal IS density: 1) Random sampling from proposal density, 2) estimating the expectation using the generated samples, and (3) updating the proposal density.

[0137] First, in the sampling step, N samples (X(t)= {^^^^^^^^}^^^^^^^=^1 are randomly drawn from the proposal density at t-th iteration step, denoted as ^^^^(^^^^) (0)^^^^. The initial proposal density, ^^^^^^^^(^^^^), is typically chosen as either the original distribution fX(x) or a uniform distribution. If prior knowledge about the region of importance is available, it can be used to construct a more informative initial proposal. For example, in seismic hazard analysis, if it is known that certain ranges of earthquake magnitude and distance contribute significantly to the hazard, the initial proposal distribution can be tailored to sample more densely within those ranges.

[0138] Secondly, in the estimation step, the integrands g(x)fX(x) in Eq. (6) are evaluated at the sampled points, and the likelihood ratios (wi) between the target density g(x)fX(x) and thecurrent proposal density q(t)X(x) are computed. These likelihood ratios serve both to estimate the expectation of g(x)and to guide the update of the proposal density for the next iteration:

[0139] The estimated value of the integral at the t-th iteration,, is obtained by taking the average of wi:Note that is unbiased as it remains an importance sampling estimator (Eq. (9)).

[0140] Third, in the updating step, the weights ^^^^(^^^^)are normalized such that ∑1. normalized represents a probability mass function proportional to the optimal density. Then, a new sample setis created by resampling from the original sample set X(t), using the weightsas the resampling probabilities. This resampling process preferentially selects samples where ^^^^ is greater than(x) to result in a set of that more closely represents the optimal sampling density ^^^^^∗ (^^^^)^^^. If the proposal density ^^^^^^^^(x) exactly matches the optimal sampling density ^^^^(^^^^)(x), all weightsbecome equal, both X(t)and Y(t)are effectively drawn from the same distribution, ^^^^∗. If the densities differ, the proposal distribution is updated by fitting the resampled points Y(t)to a presumed distribution, denoted by pX (x; η), where η is the parameters for the presumed distribution, mean and standard deviation in normal distribution, thereby obtaining an updated proposal density ^^^^(^^^^+)^^^^(^^^^).

[0141] Ideally, this iterative process continues until all weightsconverge to meaning ^^^^(^^^^)^^^^(^^^^) = ^^^^^∗^^^. However, in practice, such convergence is rarely achieved because the assumed distribution pX(x;η) may not share the same functional form as the true optimal sampling density. Instead, the iteration is typically terminated once the distributions of normalized weightsbecome sufficiently similar, indicating that updates to the proposal density are negligible. The Kolmogorov–Smirnov (K-S) 236 D statistics can be used (Kolmogorov, 1933).

[0142] After the termination of the loop, the process estimates the integral ^̂^^^(^^^^)and the interated proposal density ^^^^(^^^^)^^^^(x). In the context of PSHA,represents the mean hazard,and ^^^^(^^^^)^^^^(x) will be key to the computation of fractiles and sensitivity analysis, which will be explained below.

[0143] Before introducing the PSHA computational framework based on PMC-AIS, it is helpful to first discuss the mathematical formulation of mean and fractile hazards, and sensitivity analysis in order to establish their connection with the importance sampling framework. One should note that the individual hazard integrates the aleatory uncertainty variables, and mean and fractile hazards are obtained through the multiple runs of individual hazard calculations across the possible alternative seismic source and ground motion models, the epistemic uncertainty of the model.

[0144] By combining equations Eq. (2) and Eq. (3), the mean hazard can be expressed as:= ^^^^(^^^^)^^^^(^^^^(^^^^;^^^^) > ^^^^)^^^Φ^ |Θ(^^^^|^^^^)^^^^^^^^^^^^^^^^^^^^,^^^^

[0145] This formulation implies that in order to compute the mean hazard, it is not necessary to separate aleatory and epistemic uncertainty variables. This offers a substantial computational advantage, as it avoids the need for nested integrations—where the inner integral corresponds to individual hazard computation, and the outer integral accounts for the distribution over epistemic uncertainty. The IS formulation of mean hazard is:

[0146] The IS MC estimate of the mean hazard is:where N is the total number of samples, (Φi, Θi) are sampled from qΦ,Θ(φ, θ). The optimal sampling density which makes the variance of ^̂^^^(^^^^) zero is given by:

[0147] By rearranging Eq. (15), one obtains the following expression: (16)(54)where λ(a|θ ) denotes the hazard (i.e., exceedance rate) conditional on the epistemic uncertainty parameter θ—that is, the individual hazard given the seismic source and ground motion models. Here, λ(a) represents the mean hazard, and ^^^^Θ∗(^^^^) is the importance sampling (IS) optimal sampling density of the epistemic uncertainty variables, computed by marginalizing the IS optimal sampling density (Eq. (15)) over the aleatory uncertaintyvariables φ: ∫^^^^ ^^^^∗Φ,Θ (^^^^, ^^^^)^^^^^^^^. Here, fΘ(θ) denotes the original PDF of the epistemicuncertainty variables. By repeatedly calculating Eq. (16) using random Θ samples following fΘ(θ), a list of individual hazards can be obtained.

[0148] Recall that the CDF of the hazard, from which fractile hazards can be derived, must be computed empirically by sorting the individual hazards, λ(a|θ ), computed from numerical integration of Eq. (1) (see Eq. (4) and Eq. (5)). The Eq. (16) eliminates the need for numerical integration typically required for individual hazard calculation. That is, individual hazard can be efficiently calculated if the mean hazard, λ(a), and the optimal sampling density of the epistemic uncertainty variables, ^^^^Θ∗(^^^^), are known.

[0149] Equation Eq. (16) also yields an important insight regarding the optimal sampling density and the sensitivity of each epistemic uncertainty variable to the hazard. Suppose one is interested in assessing the sensitivity of the i-th epistemic variable, θi, to the hazard. If variations in θihave a negligible effect on the hazard, then the individual hazard remains approximately constant across different values of θi, i.e., λ(a|θi) ≈ λ(a) for all θi. In this case, the marginalized optimal sampling density(^^^^^^^^)should closely resemble the original distribution Conversely, if ^^^^^^^^exerts a substantial influence on the hazard, then ^^^^∗(θ) will significantly deviate fromTherefore, the degree of deviation between marginalized optimal and original densities can serve as an indicator of the sensitivity of the hazard to the corresponding epistemic uncertainty variable. The hazards with respect to the various values of ith epistemic uncertainty variable can be calculated by reformulating Eq. (16) in terms of the i-th variable θi:where and ^^^^∗Θ^^^^ (^^^^^^^^) =∫^^^^^^^^ ^^^^∗Φ,Θ (^^^^,^^^^)^^^^^^^^^^^^^^�^^^^^^ denotes the set of epistemic uncertaintyvariables excluding�^^^^^^^^.

[0150] The variance contribution of the i-th epistemic uncertainty variable to the total hazard variance, denoted Ci, can be computed as:where ^^^^^^^^^^^^[^^^^^^^^(^^^^|^^^^^^^^]is the variance of the hazard when θi is varied, and ^^^^^^^^^^^^[^^^^(^^^^|^^^^]is thevariance of the hazard when all the epistemic uncertainty variables (θ) are varied. Although the computation of λi(a|θi) also inherently involves numerical integration, the use of Eq. (17) significantly accelerates this process. Interestingly, Ci corresponds to the first-order Sobolindex of θi, under the assumption that θi is independent of�^^^^^^^^—i.e., θi is propagated throughall other variables. The first-order Sobol index quantifies the direct contribution of each variable to the total variance of the estimate and ranges between 0 and 1. A larger value of Siindicates that the hazard is more sensitive to θi. The sum of Ci, denoted^^^^^^^^, where dθ is the number of epistemic uncertainty variables, does not necessarily equal 1. The residual term1 − ∑^^^^^^^^^^^^ ^^^^^^^^Ci represents the interaction effects among the variables.

[0151] While Eq. (16)-Eq. (18) are useful in fractile hazard calculation and sensitivity analysis. However, there are three primary challenges in utilizing the equations. First, the mean hazard λ(a) must be known. In the traditional PSHA framework, individual hazards are first computed and then averaged to obtain the mean hazard. In contrast, equation Eq. (16) requires the mean hazard to compute the individual hazard, which may seem paradoxical.Second, ^^^^Θ∗ (^^^^)) (^^^^Θ∗^^^^(^^^^^^^^))) can be derived through marginalization of ^^^^Φ∗,Θ(^^^^,^^^^)and theaccuracy of obtained fractile is directly related to how we can closely approximate the optimal sampling density. However, obtaining ^^^^Φ∗,Θ(^^^^,^^^^)is non-trivial problem, requiring well-designed algorithm to approximate this. Third, even if the optimal sampling density,^^^^Φ∗,Θ(^^^^,^^^^), becomes available, obtaining its marginal distribution, ^^^^Θ∗ (^^^^)or ^^^^Θ∗^^^^(^^^^^^^^), stillrequires numerical integration, which is computationally expensive in general. In the following section, we propose a framework to address these challenges using Gaussian Population Monte Carlo Adaptive Importance Sampling (G-PMC AIS).

[0152] PMC AIS is a general framework for numerically computing the mean and identifying a proposal IS density that closely approximates the optimal sampling density. For a target ground motion, a, applying PMC AIS to PSHA iteratively yields both the mean hazard, λ(a),and the iterated proposal IS density, ^�^^^Φ∗,Θ, which approximates the optimal sampling density ^^^^Φ∗,Θ(^^^^,^^^^). This resolves the first two challenges associated with applying Eq. (15) to obtain individual hazard estimates, as noted in the previous paragraph. Furthermore, by employing a joint, or multivariate, Gaussian distribution as the assumed distribution in PMC AIS, meaning pX (x; η) = πX(x;µ,Σ), where µ and Σ are the mean vector and covariance matrix, respectively. This addresses the third challenge—computational difficulty in marginalization. For a joint Gaussian, marginalization is straightforward: one can simply select the relevant entries from the mean vector and covariance matrix corresponding to the variables of interest, without requiring numerical integration. Additionally, the use of Gaussian distributions provides further advantages: they are easy to sample from, their PDFs are straightforward to evaluate for weight calculations, and they effectively capture dependencies among variables.

[0153] The embodiments here refer to this framework as Gaussian Population Monte Carlo Adaptive Importance Sampling PSHA (G-PMC AIS PSHA). The framework comprises three primary stages: 1) Initialization, 2) Estimation of the mean hazard and the iterated proposal IS density, and 3) Fractile calculation, shown in the Algorithm 2, and FIG.19).

[0154] Stage 1. Initialization In Line 1 of Algorithm 1, the initialization stage begins by setting the iteration step t = 0 and initializing the Kolmogorov–Smirnov (K-S) D statistic d = 1, the maximum possible value, to ensure at least one iteration. The stopping criterion is defined by ε = 0.1, explained below. The K-S D statistic measures whether two probability distributions are likely to come from the same underlying distributions.

[0155] Next, the joint probability density of φ, aleatory variables such as magnitude, distance, and ground motion residual, and θ, epistemic variables such as the b-value, mmax, and alternative ground motion models, is defined (Line 2). The algorithm defines it as theconditional density of φ given θ as^^^^). This process is equivalent to constructing thelogic-tree structure and defining each end branch’s PDF of magnitude, distance, and ground motion distributions in conventional PSHA studies. Note that, unlike traditional approaches, the distribution of epistemic uncertainty variables in this framework can be continuous and does not require discretization. A simple two-dimensional example involving one aleatory variable (φ: magnitude M) and one epistemic variable (θ : ∆µ, median ground motion model) is depicted in the "Stage 1. Initialization" box in FIG.19.

[0156] Finally, the initial proposal density ^^^^(0)Φ,Θ(^^^^,^^^^)is selected (Line 3). Assuming a joint Gaussian form, the mean vector and covariance matrix must be specified. If prior knowledgeexists regarding important regions of each variable for the target ground motion a, the joint Gaussian can be centered on those values. Otherwise, a common strategy is to center the mean at the midpoint of the integration domain and use a covariance matrix with large variances and zero covariances to span the full range. These parameters will be updated in subsequent iterations to reflect the important regions and interdependencies.

[0157] Stage 2. In this stage, the process estimates both the mean hazard ^�^^^(^^^^) and iteratedproposal IS density, which approximates the optimal sampling density. This is an iterative process, and considers the t-th iteration step.

[0158] A random sample set ^⃗^^^(^^^^)consisting of φ and θ is drawn from. The indicator function I(u)(^⃗^^^(^^^^)) > a) is evaluated for each sample, and the results are multiplied by the earthquake occurrence rate ν(θ) and the likelihood ratio (importance weight) to obtai ^^^^n the weight vector ^^^⃗^ (Lines 5–6). The mean hazard for t-thiteration, ^�^�^^⃗ (^^^^)(^^^^), is then estimated as the average of ^^^⃗^ ^^^^ (Eq. (14), Line 7).

[0159] Next, the process updates the proposal density to better approximate the optimal sampling density. A new sample set ^�^^⃗^ (^^^^) is formed by resampling from ^⃗^^^(^^^^) according to^^^^ is proportional to the optimal sampling density,by the proposal density in the current iteration. The res(^^^^)ampling process redistributes the original samples ^⃗^^^ to follow optimal sampling density. The updated proposal density ^^^^(^^^^+1)Φ,Θ (^^^^, ^^^^). is then obtained byfitting a joint Gaussian distribution to ^�^^⃗^ (^^^^) via maximum likelihood estimation (Lines 9–10).

[0160] To decide whether to terminate or continue the iteration, the K-S D statistic is evaluated (Line 11). This statistic measures the maximum difference between the CDFs ofA smaller d indicates convergence, and iteration terminates when d < ε, where ε = 0.1 in these embodiments. The outputs of this stage are the mean hazard estimate ^̂^^^(^^^^) and the approximated optimal (iterated proposal) density, ^^^^) (Line 13). An example is shown in the "Stage 2" box of FIG. 19,the evolved proposal density captures the integration’s important region near φ ≈ 5.7 and θ ≈ 0.2.

[0161] Stage 3. In this stage, ∗^^^^) is marginalized to obtain, ^�^^^Φ,Θ (Line 14). Thisdone by extracting the θ-related elements of the mean vector and covariance matrix without additional computation. Random samples of θ are drawn from the original epistemicuncertainty distribution fΘ(θ) to form the sample set ^^^^ (Line 15). These are visualized in the leftmost plot of the "Stage 3" box in FIG.19, where the solid curve represents fΘ(θ) and dots indicate the samples representing key fractiles (1st, 16th, 50th, 84th, and 99th percentiles).

[0162] Next, likelihood ratios between ^�^^^Φ∗,Θ(^^^^,^^^^) and fΘ(θ) are computed, shown in the second plot in FIG.19. These ratios are scaled using the mean hazard ^̂^^^(^^^^), producing individual hazard estimates ^^^^(a|θ) (Eq. (16), Line 16). The pth fractile hazard, ^̂^^^^^^^(^^^^), is then identified such that p% of the individual hazards are less than or equal to λp(a) (Eq. (5), Line 17). This process is illustrated in the final plots of the "Stage 3" box in FIG.19. This overall procedure (Stages 1 to 3) is repeated for each target ground motion intensity of interest.

[0163] It is important to note that the overall computational effort is dominated by “Stage 2” due to repeated evaluations of the marginal hazard for numerical integration. In contrast, “Stage 3”, fractile calculation, is computationally lightweight, involving only likelihood ratio computations between known distributions. Table 4:

[0164] The discussion now considers two different seismic source configurations: areal and fault sources, shown in FIGs.20A-20B. These are commonly used seismic source types in practical PSHA applications. The earthquake occurrences are modeled using a truncated exponential distribution, with a minimum magnitude mmin set to 5.0. The annual occurrence rate for earthquakes with magnitudes exceeding mmin, denoted as ν, is set to 0.01. The ground motion model is employed that assumes a VS30of 760 m / s.

[0165] This process incorporates four epistemic uncertainty variables: two related to the seismic source—the b-value and mmax—and two associated with the ground motion model— its median and standard deviation. The b-value and mmaxare assumed to follow normaldistributions: b ∼ ^^^^(1.0, 0.1) and mmax ∼ ^^^^ (7.0, 0.3), with truncation to the ranges [0.7,1.1] and [5.9, 7.1], respectively, to prevent unrealistic values. Uncertainty in the naturallogarithm of the median ground motion is represented as ∆µ ∼ ^^^^ (0, 0.1), corresponding toapproximately a 10% variation in peak ground acceleration (PGA) within one standard deviation. Similarly, the uncertainty in the standard deviation of the natural logarithmicground motion is modeled as ∆σ ∼ ^^^^ (0, 0.05).Algorithm 2: G-PMC AIS PSHA computational framework pseudocode given target ground motion a [1. Initialize] 1: t = 0, d = 1, ε =0 / 1 2: Defneprobability of se(^^^^)lecting ^^^^^^^^is12: ^^^^ ← ^^^^ + 1[3. Fractile] 14:15: ^^^^ = {Θ^^^^}^^^^^^^^=1 , ^^^^~ ^^^Θ^ (^^^^)16:17:18: return

[0166] The process estimates hazard levels at 0.13g, 0.32g, 0.64g, and 1.1g, which correspond to annual exceedance probabilities of approximately 10−3, 10−4, 10−5, and 10−6for the fault source, and ~ 1.3 × 10−4, 1.2 × 10−5, 1.2 × 10−6, and 1.3 ×10−7for the areal source.

[0167] The discussion now turns to comparing the proposed framework, “G-PMC AIS,” with logic-tree approaches using widely adopted discretization schemes: three- and five-point approximations, denoted as “LT(3)” and “LT(5)”, respectively in Table 3. Given four epistemic uncertainty variables, the total number of end-branches (Nθ) is 81 and 625 for the “LT(3)” and “LT(5)”, respectively. The total number of marginal hazard evaluations is Nt = Nθ × Nφ , where Nφ is the number of samples used per individual hazard computation. For both methods, individual hazard estimation is performed using the framework above, which is significantly faster than traditional techniques. The process uses Nφ = 104, which ensureshigh accuracy of estimates (COV ≤ 2.5%) for most cases considered.

[0168] Additionally, the performance of a continuous distribution based conventional Monte Carlo (“MC”) method is evaluated. This approach is similar in structure to the logic-tree method but draws samples directly from the continuous distributions without approximation. For consistency, the same Nφ= 104 is used. It guarantees the true fractile if sufficient number of individual hazards are estimated. The benchmark CDF is derived from 1,000 MC realizations.

[0169] To evaluate the accuracy of mean hazard estimates, COV is used, representing the relative variability of the estimate with respect to mean. A lower COV indicates a more robust estimate. For fractile hazard accuracy, the Kolmogorov–Smirnov (K-S) D 384 statistic was used to compare the CDF of each method’s fractile hazard with a benchmark CDF derived from 1,000 MC realizations. Note that a K-S D statistic below 5% generally indicates good agreement. For instance, if hazard is assumed log-normally distributed, a 5% K-S D corresponds to a mean difference of only 1%–2.3% depending on the assumed standarddeviation. Computational efficiency is measured by the total number of marginal hazard computations, Nt.

[0170] FIGs.21A-21B shows the coefficient of variation (COV) of mean hazard estimates as a function of Nt, and FIGs.22A-22B presents the Kolmogorov–Smirnov (K-S) D statistics for the fractile hazards. Across all scenarios—both areal and fault sources and for all ground motion levels—the G-PMC AIS framework demonstrates superior efficiency.

[0171] To achieve a 1% COV for PGAs of 0.13g, 0.32g, 0.64g, and 1.1g, G-PMC AIS requires 29, 29, 21, and 16 times fewer samples than the “LT(3)”, and 224, 224, 162, and 126 times fewer than “LT(5)”, respectively, shown in FIGs.21A-21B. Compared to “MC”, the sample reductions are even more pronounced: 286, 559, 1666, and 3,775 times fewer samples are required by G-PMC AIS, respectively. This underscores the inefficiency of MC method at higher ground motion levels, where rare-event modeling becomes computationally burdensome, unlike G-PMC AIS, which maintains stable performance.

[0172] Although G-PMC AIS excels in mean hazard estimation, high accuracy in the mean does not inherently guarantee accurate fractile estimates. To evaluate this, K-S D statistics were computed at Nts corresponding to 1% COV in mean hazard estimation, 400 achieved using G-PMC AIS. The resulting K-S Ds for 0.13g, 0.32g, 0.64g, and 1.1g are 5.6%, 4.0%, 2.8%, and 3.0%, respectively, 401 indicating good agreement with the benchmark CDFs shown in FIGs.22A-22B.

[0173] While the five-point logic tree, LT(5), yields more accurate fractile estimates than the three-point, LT(3), as seen in FIGs.18A-18C, its computational cost increases significantly with the number of epistemic variables. Thus, although preferred for accuracy, it suffers from a computational burden due to the rapid growth in Nt.

[0174] An interesting contrast is observed in the performance of the MC method between mean and fractile estimates, FIGs.21A-21B and 22A-22B. For the mean, performance deteriorates with increasing ground motion level due to the rarity of higher threshold exceeding events. However, in the case of fractiles, the K-S D statistic is relatively insensitive to ground motion, since all samples contribute equally to the CDF regardless of event rarity.

[0175] It is important to note that the efficiency of LT(3), LT(5), and MC methods depends heavily on the number of samples used or individual hazard calculations. While the processes use Nφ= 104, reducing this to 103—for example, corresponding to approximately 20 magnitude bins (e.g., increments of 0.1 from 5 to 7) and 50 spatial samples—still results inG-PMC AIS outperforming all other methods. Moreover, optimizing the grid size for such methods demands significant prior knowledge or sensitivity analysis, incurring additional computational cost. In contrast, G-PMC AIS operates effectively without such prerequisites.

[0176] For fault sources, the performance gap remains substantial. At a 1% COV, G-PMC AIS is 13–24 times, 98–183 times, and 123–2268 times more efficient than the three-point logic tree, five-point logic tree, and MC methods, respectively. The K-S D 418 values for the 1% COV cases are 3.6%, 2.8%, 3.9%, and 2.8% for 0.13g, 0.32g, 0.64g, and 1.1g, respectively, again demonstrating strong agreement with the benchmark. Interestingly, when using small sample sizes, fault sources outperform areal sources in 420 fractile estimation, likely due to reduced problem dimensionality, shown in FIGs.22A-22B. In the implementation of the embodiments, one spatial variable (longitude) is eliminated in fault source, due to the vertically lying fault, thereby reducing computational complexity.

[0177] These numerical results demonstrate that the proposed G-PMC AIS framework is, in general, faster than any of the alternative methods considered. Accurate fractile hazard estimation depends on whether the iterated density adequately approximates the optimal sampling density. In practical PSHA applications, this approximation quality is not known, as verifying each point estimate against the true value would require the same computational effort as traditional numerical integration, negating the benefit of the proposed method. Nevertheless, the numerical examples suggest that when a sufficiently low COV (e.g., < 3%) is achieved for the mean hazard, the corresponding fractiles also exhibit good agreement.

[0178] Conversely, when the number of samples is insufficient (e.g., N = 103 in areal source case; FIG.21A, the K-S D statistic exceed 50% in some cases, indicating a significant deviation from the true fractile and illustrating a failure to approximate the optimal sampling density under limited sampling. This highlights the importance of using an adequate number of samples in G-PMC AIS, though the required sample sizes remain substantially lower than for other methods. Unlike MC or logic-tree approaches, where accuracy improves incrementally with sample size, G-PMC AIS exhibits a “phase transition” behavior— convergence toward the optimal sampling density occurs only after a critical threshold in Nt is surpassed.

[0179] For the areal source case, a sensitivity analysis based on the likelihood ratio between the marginalized densities, as defined in Eq. (18), is presented in FIG.23. These results closely match the benchmark first-order Sobol indices with a maximum difference of 0.013. The results indicate that uncertainties in the ground motion model—specifically, the median(∆µ) and standard deviation (∆σ)—are the dominant contributors to overall sensitivity, consistent with findings from previous PSHA studies. In contrast, the influence of seismic source parameters, such as mmaxand the b-value, is relatively minor. The notable shift in the distribution of ∆σ—highlighting its central role in shaping the hazard—can be seen in the comparison between the marginalized iterated proposal densities (approximated optimal sampling density) and the original prior distributions.

[0180] Additionally, the interaction effect among variables, represented by ρ in FIG.23, generally exhibits minimal contribution. Probabilistic Seismic Hazard Analysis (PSHA) has traditionally been a computationally intensive task, particularly when complex logic-tree structures with numerous epistemic uncertainty variables are involved.

[0181] The embodiments demonstrated that approximating the underlying continuous distributions with discrete representations can inherently lead to biased estimates of both mean and fractile hazards, especially when arbitrary weighting schemes are used in the approximation process. From a simple numerical example, it is shown that the mean, 16th,and 84th fractile hazards differ by factors of ∼ 1.8, 2.7, and 1.8, respectively, which can beeven larger when the lower exceedance probability is considered. Furthermore, the computational burden of hazard calculation increases exponentially with the number of epistemic uncertainty variables.

[0182] To address these issues, the embodiments may use continuous distributions to preserve the original characteristics of the epistemic variables. Furthermore, to mitigate the computational challenges associated by incorporating high-dimensional continuous distributions, the embodiments propose a PSHA computational framework based on Gaussian Population Monte Carlo Adaptive Importance Sampling (G-PMC 456 AIS).

[0183] The key insights of the G-PMC AIS approach in calculating PSHA are: 1) the mean hazard can be efficiently estimated G-PMC AIS approach based on the formulation treating aleatory and epistemic uncertainty variables jointly, thereby avoiding the need for double- nested integration and improving computational efficiency; 2) individual hazard estimates used for fractile computation can be obtained efficiently, almost without additional cost, by utilizing the likelihood ratio between marginalized optimal sampling density, which is a by- product of mean estimation, and the original density, with the easy marginalization by adopting joint Gaussian distribution; 3) The hazard sensitivity analysis also can be done easily in the same way.

[0184] The embodiments demonstrated the effectiveness and validity of the proposed methodology through numerical examples. The results indicate that the G-PMC AIS framework outperforms three-point approximation logic tree, five-point approximation logic tree, and continuous distribution based MC approach by a factor of up to 29, 224, and 3775 times while achieving high accuracy (COV < 1%) of mean hazard estimation.

[0185] Moreover, K-S D statistic between fractile hazards for those cases with high meanhazard accuracies and the benchmark all fall below ∼ 5%, indicating extremely similarhazard fractile distribution, indicating the suggested framework is able to approximate 469 the optimal sampling density well. The likelihood ratio between marginalized each epistemic uncertainty variable’s optimal sampling density and marginalized presumed epistemic uncertainty also enables us to quantify the sensitivity of each variable to the hazard. The embodiments can be effectively applied to a wide range of PSHA applications involving complex logic-tree structures with extensive epistemic uncertainties.

[0186] All features disclosed in the specification, including the claims, abstract, and drawings, and all the steps in any method or process disclosed, may be combined in any combination, except combinations where at least some of such features and / or steps are mutually exclusive. Each feature disclosed in the specification, including the claims, abstract, and drawings, can be replaced by alternative features serving the same, equivalent, or similar purpose, unless expressly stated otherwise.

[0187] Additionally, this written description makes reference to particular features. It is to be understood that the disclosure in this specification includes all possible combinations of those particular features. For example, where a particular feature is disclosed in the context of a particular aspect, that feature can also be used, to the extent possible, in the context of other aspects.

[0188] Also, when reference is made in this application to a method having two or more defined steps or operations, the defined steps or operations can be carried out in any order or simultaneously, unless the context excludes those possibilities.

[0189] Although specific aspects of this disclosure have been illustrated and described for purposes of illustration, it will be understood that various modifications may be made without departing from the spirit and scope of the invention. Accordingly, the invention should not be limited except as by the appended claims.

Claims

WHAT IS CLAIMED IS:

1. A computer-implemented method of performing probabilistic seismic hazard analysis (PSHA), comprising: receiving a seismic source model that provides a probability distribution of magnitude and distance from a site of interest to where earthquakes may occur, and a ground motion model that provides a probability distribution of an intensity of ground motion that will be observed at the site of interest; constructing a joint probability distribution from the seismic source model and the ground motion model; using non-parametric adaptive importance sampling of the joint probability distribution to obtain an optimal sampling density; computing a PSHA hazard estimate for the site of interest; computing a hazard deaggregation distribution for the site of interest; and providing the PSHA hazard estimate and the hazard deaggregation to a user.

2. The computer-implemented method as claimed in claim 1, wherein using the non- parametric adaptive importance sampling comprises: defining an initial sampling density as a uniform sampling density; dividing the data into grids and assigning a probability to each grid; iteratively adjusting the probability density for each grid depending upon a contribution by that grid to a hazard estimator until no further improvement of the hazard estimator occurs.

3. The computer-implemented method as claimed in claim 2, wherein the grids comprise cuboids for three-dimensional data.

4. The computer-implemented method as claimed in claim 1, wherein the seismic source model comprises a single source model.

5. The computer-implemented method as claimed in claim 1, wherein the seismic source comprises a combined source model.

6. The computer-implemented method as claimed in claim 5, wherein a source probability for each source in the combined source model is implemented in the joint probability density function.

7. The computer-implemented method as claimed in claim 5, wherein the PSHA mean hazard estimate is computed for each source, and the estimate for each source is summed with the other sources.

8. The computer-implemented method as claimed in claim 1, wherein computing a hazard deaggregation curve comprises using the optimal sampling density as deaggregation density.

9. A computer-implemented method of performing probabilistic seismic hazard analysis (PSHA), comprising: receiving a seismic source model that provides a probability distribution of magnitude and distance from a site of interest to where earthquakes may occur, and a ground motion model that provides a probability distribution of an intensity of ground motion that will be observed at the site of interest; constructing a joint probability distribution of epistemic and aleatory uncertainty variables from the seismic source model and the ground motion model; estimating a mean hazard for each of a series of iterations by: drawing a random sample set from a proposal sampling density; evaluating an indicator function for each sample and multiplied by an earthquake occurrence rate and a likelihood ratio to produce a weight vector. estimating a mean hazard from an average of the weight vector;updating the proposal sampling density to a next iteration proposal sampling density; and terminating the series of iterations when the proposal sampling density and the next iteration proposal sampling density converge, the series of iterations producing a mean hazard estimate and the next iteration proposal sampling density as an approximated optimal sampling density.

10. The computer-implemented method as claimed in claim 9, further comprising: extracting random epistemic uncertainty variables from the approximated optimal sampling density; applying a likelihood ratio to each random sample; scaling each likelihood ratio by the mean hazard estimate to produce individual hazard estimates; and sorting the individual hazard estimates to produce a hazard fractile (or percentile).

11. The computer-implemented method as claimed in claim 10, wherein the epistemic uncertainty distribution is continuous.

12. The computer implemented method as claimed in claim 9, wherein updating the proposal sampling density comprises: resampling the random sample set according to weights in the weight vector to produce a new sample set; fitting a joint Gaussian distribution to the new sample set to produce the next iteration proposal sampling density.

Images