Patent ID: 12236369

DETAILED DESCRIPTION

As described in various embodiments, a system, method, and corresponding non-transitory computer readable medium storing machine interpretable instructions (e.g., articles of manufacture directed to computer program products) are directed to an improved computational approach to adaptive importance sampling for Monte Carlo (MC) estimation on high dimensional integrals.

MC methods are a class of powerful algorithmic computational approaches that allow estimation of high-dimensional integrals or sampling from nontrivial probability distributions. MC estimations can employ repeated random sampling of target distributions to determine parameters associated with the target distribution. For example, MC estimation may be used to estimate an expected value of a probability distribution.

Experimentally, traditional approaches in conducting MC analyses were found to be insufficient for certain use cases, especially where tail behavior in the target distribution required accurate capture.

Applicant proposes a method of conducting MC analyses based on adaptive importance sampling (AIS). AIS includes two separate features; (1) importance sampling, wherein a target distribution is modelled with an importance distribution, wherein the importance distribution overemphasizes, or is biased, sampling values from an area of the target distribution which is considered to be more important (e.g., sampling to a greater extent from tail regions of a distribution indicative of financial solvency, where the tail represents insolvency) and (2) adaptive importance sampling, where the importance distribution is adjusted over time to more closely match the target distribution based on samples generated from the initial importance distribution.

To generate the MC estimation, adaptive importance sampling (hereinafter referred to generally as importance sampling) draws samples from the latest importance distribution, which samples are used to determine a weight (e.g., a ratio) of the target distribution probability density at the sample and the importance distribution probability density at the sample (i.e., the importance weight). The importance weight is subsequently used to construct a MC estimate.

The following examples are presented here and for reference and to provide an illustration of the uses of adaptive importance sampling and MC estimation.

In an example embodiment, the proposed system generates MC estimates representative of a likelihood that a hospital or other medical service establishments will have a certain amount of infected staff at a particular time. The importance distribution (alternatively referred to as the proposal distribution) may be used to represent one or more of the following likelihoods (the target distribution); the likelihood an infected patient infects a staff member, the likelihood that a patient is infected, the likelihood that a staff member is infected, the likelihood that an infected staff member infects a patient, and so forth. The aforementioned target distributions may be the subject of pre-existing studies, such that at least some characteristics of the target distribution of said phenomena in previous circumstances are known or suspected based on experimental data (e.g., the likelihood that a patient infects a staff member may be suspected to be represented by a ten dimensional Bi-modes diagonal Gaussian distribution based on previous data of infectious diseases).

The proposed model may, via adaptive importance sampling, continually adapt a proposal distribution (e.g., a Student's T distribution) to more accurately represent an estimated ten dimensional Bi-modes diagonal Gaussian distribution. The proposed adaptive importance sampling adapts the proposal distribution by generating samples from the proposal distribution, and assessing them relative to the previous experimental data.

The generated samples may subsequently be used to determine importance weights (as described herein) which are used as samples in a MC estimation.

In another example embodiment, the proposed system generates MC estimates representative of a likelihood of production line failure. Importance distributions may be used to represent the following target distributions: the likelihood of a particular piece of equipment malfunction, the likelihood of staff absences, the likelihood of material shortages, and so forth. The likelihood of the particular piece of equipment may be suspected based on previous information provided by a manufacturer of said equipment, the likelihood of staff absences may be derived from existing staff absence data, and so forth.

As with the hospital example, the proposed system models the target distributions through adaptive importance sampling, and uses at least some samples drawn from the continually updated importance distributions to generate a MC estimate.

In a further example embodiment, the proposed system generates MC estimates representative of a likely price of a particular financial instrument. Importance distributions may be used to represent the following likelihoods; the likelihood of consumer sentiment shifting, the likelihood and magnitude of interest rates changing, and so forth.

An approach to generating MC estimates is Markov Chain MC (MCMC), where samples from a target distribution are obtained by construction and running of a Markov chain for a long time. However, the MCMC framework suffers from a number of technical limitations. First, it is an inherently serial construction, limiting the ability of MCMC methods to exploit the parallel structure of modern computational hardware (e.g., GPUs which are designed to perform homogenous computations in parallel). Second, samples generated with the MCMC method that are temporally close are highly correlated with each other, either increasing the amount of time the model needs to operate to produce an accurate MC estimation, or underrepresenting less explored portions of the target distribution in the final MC estimation. Finally, the constraints placed on Markov Chain transition kernel to implement MCMC make them difficult and often unintuitive to design.

The proposed approach overcomes the serial nature of MCMC by implementing importance sampling to take advantage of parallelizable computing infrastructure. Samples are independently generated (e.g., samples are generated without reference to a previous sample) and evaluated. Furthermore, as a result, the approach is unbiased and guarantees asymptotic convergence.

A technical challenge associated with implementing importance sampling is that the choice of importance distribution may be incorrect, or a poor approximation of the target distribution, such that the variance of the MC estimator can be arbitrarily high or infinite. As a result of these technical challenges, many approaches discard importance sampling as a framework for general purpose MC inference (e.g., arguing that if enough is known about a target distribution to construct a good importance distribution, then there is likely little left to be done).

A further technical challenge associated with implementing importance sampling is that the importance distribution should be heavier tailed than the target distribution. Higher order moments are sensitive to the tail behaviour of the distribution which can lead to high variance importance sampling estimators if that behaviour is not well captured.

Some approaches construct methods that automatically adapt mapping of an importance distribution to the target distribution. The methods, called adaptive importance sampling have been limited by the simplistic nature of their importance distributions which are adapted. For example, they have largely been constructed either as a single Gaussian or as mixtures of Gaussians.

Applicant proposes an approach whereby normalizing flows, a new class of flexible parametric probability distributions, are used to overcome the technical challenges associated with both selecting and generating an importance distribution heavier tailed than the target distribution, and generating complex importance distributions for AIS. Normalizing flows are a framework for constructing complex probability distributions by applying a succession of bijective transformations to a random variable with a simple, known distribution. Normalizing flows represent a potentially powerful class of models for density estimation. The use of normalizing flows in solving some of the technical problems described herein is new, as the study of normalizing flows is in a nascent phase.

The proposed approach incorporates normalized flows to further configure adaptive importance sampling to issues that arise in more general MC estimation for statistical inference, and not just to limited cases. Alternatively stated, the proposed approach constructs a framework of adaptive importance sampling to perform MC estimation on challenging target distributions (e.g., target distributions with heavy tails) using normalizing flows to represent the importance distribution.

As described in some embodiments, a novel loss function is implemented that specifically aims to minimize the variance of the importance sampling estimator for higher order moments. This encourages the machine learning architecture to explore and capture the tail of the target distribution while still learning the general shape and typical set.

The machine learning architecture's ability to learn heavy tailed distributions is explored empirically. Such distributions can be used in respect of in the Hierarchical Bayesian model. In Bayesian inference, a fairly weakly informative prior distribution is desired.

Unfortunately, the empirical results show that some popular normalizing flow models with standard Gaussian as base distribution struggle to model the heavy tailed target distributions. To address this technical challenge, Applicants, in some embodiments, propose an approach to replace the Gaussian base distribution in the normalizing flow with a Student-t distribution. In example embodiments, applicant proposes a base or initial importance distribution to address the variance issue such as heavy tailed distributions including Half-t distribution, Half-Cauchy and log-Normal. Such choice of base or initial importance distributions generally leads to a heavy tailed resulting importance distributions. Moreover, the importance distribution can be a generalized Gaussian distribution which can be parameterized. This provides additional parameters, the degree of freedom, which are learned along with the other parameters of the flow. Further, it also allows the estimated distribution to better capture heavy tailed distributions.

Referring now toFIG.1, a diagram of a system for computing MC estimates is disclosed. Computing device100of the system includes processor102, databases122, and machine learning architecture110.

In example embodiments, the computing device100is connected, via networks116, to external databases118or external devices120to retrieve information related to generating the MC estimate. For example, computing device100may be configured to retrieve samples, distribution models, or other model inputs for either the target distribution or the importance distribution (as described herein) from external databases118or external devices120. According to some embodiments, for example, the computing device100is in communication with database122of computing device100to retrieve samples, distribution models, or other model inputs.

The processor102is in communication machine learning architecture110, including sample generator106, and MC estimator104, to generate samples based on the importance distribution, and to generate the MC estimate based on the generated samples. In example embodiments, sample generator106and MC estimator104are non-transitory computer readable machine code stored in memory of computing device100, which when executed by the processor cause a computing device100to perform any of the methods described herein.

The machine learning architecture110includes a plurality of parameters which are used to define the sample generator106and MC estimator104. For example, the sample generator may be a function defined by one or more parameters whose values are iteratively updated with gradient based methods and a loss function. The selection of loss function defines in part which of and how the value of the parameters change with update iterations. For example, some update iterations will only update a few parameters a small amount, and other iterations may update all parameters a large amount. The values of the parameters of machine learning architecture110therefore define in part the MC estimates generated by computing device100.

MC estimator104generates MC estimates based on one or more sample observations (hereinafter referred to as samples for ease of reference) obtained by the sample generator106, for example via estimate generator114. MC estimator104also determines whether samples generated by the sample generator106belong to the target distribution with error estimator112. In example embodiments, as described herein, the MC estimator104operates according to different states, wherein in a first state the MC estimator104exclusively evaluates error, and in a second state the MC estimator104both evaluates error and accumulates samples to generate the MC estimate.

Estimate generator114repeatedly, and potentially in parallel, evaluates random samples obtained from the sample generator106to obtain numerical results. For example, some example embodiments of estimate generator114of MC estimator104are described below:

Let {tilde over (p)}(x) be a non-negative, real-valued function defined inDsuch that ∫D{tilde over (p)}(x)dx=Z is a finite real positive value.

Then set {tilde over (p)}(x) to be an unnormalied probability distribution (i.e., the target distribution) such that

p⁡(x)=1Z⁢p˜(x)
is its normalized form and Z is called the normalizing constant of {tilde over (p)}(x).

Estimate generator114, given {tilde over (p)}(x) or p(x), would like to estimatep(x)[ƒ (x)].

Estimate generator114, uses the samples generated by the sample generator106based on the importance distribution (denoted by q(x(i)) or simply q(x)) as samples of the target distribution and thereafter evaluates the samples to determine a weight (w(i)) and aggregates the weights to determine an expected value Σiw(i)ƒ(x(i)), or other parameter, of the target distribution p(x) based on the following relations:

𝔼p⁡(x)[f⁡(x)]=∫f⁡(x)⁢p˜(x)Z⁢d⁢x=∫f⁡(x)⁢p˜(x)Z⁢q⁡(x)q⁡(x)⁢dx=𝔼q⁡(x)[f⁡(x)⁢p˜(x)Z⁢q⁡(x)]≈1N⁢∑iN⁢f⁡(x)⁢p˜(x(i))Zˆ⁢q⁡(x(i))⁢where⁢(Zˆ=1N⁢∑iN⁢p˜(x(i))q⁡(x(i))).Let⁢W(i)=p˜(x(i))q⁡(x(i))⁢and⁢w¯(i)=w(i)∑i⁢w(i).Then⁢⁢Zˆ=∑i⁢w(i)N⁢and⁢𝔼p⁡(x)[f⁡(x)]≈∑i⁢w¯(i)⁢f⁡(x(i)).

The samples are thereafter aggregated based on (with reference toFIG.1B).

For a sequence of importance distributions {q(x; θt)}t=0Twhere T is the total number of iterations and q(x; θt) is the importance distribution parameterized by θtat iteration t.

At block1002, for each q(x; θt), the approach includes drawing a set of samples {xt(i)}˜q(x; θt), where θtis determined by samples from the previous iteration {xt-1(i)}.

At block1004, the weights are determined and normalized via:

w¯t(i)=wt(i)∑t,i⁢wt(i).

At block1010, a MC estimation is performed based on:p(x)[ƒ(x)]≈Σt,iwt(i)ƒ(xt(i)).

Once the samples have been aggregated, or prior to the samples being aggregated to determine an expected value (e.g., where the MC estimator104operates according to a first stage) the MC estimator104may determine an error of the generated sample with error estimator112. For example, as will be described below, the error estimator112may employ a loss function to determine an error associated with the samples.

Error estimator112of MC estimator104may implement adaptive importance sampling, where the importance distribution of importance distribution modeller108A used to obtain the samples is updated based on previous samples. In example embodiments, error estimator112may simply determine the weights of the generated samples and a separate gradient determiner (not shown), or other component may implement adaptive importance sampling based on the values determined by error estimator112.

In order to update the plurality of parameters of machine learning architecture110of the sample generator106(e.g., the plurality of parameters that define the importance distribution modeller108A), the machine learning architecture110is configured with a loss function defined to processes samples previously generated by the sample generator106. The plurality of parameters are updated based on a gradient of the loss function.

In previous approaches, the error estimator112would update the parameters θtof the mapping function (e.g., at least some of the plurality of parameters), with a forward KL divergence as a loss function to determine the extent of the updates to the parameters, which frame work is shown below:

OK⁢L(θ)=KL[p˜(·)||q⁡(·;θ)](1)=∫p˜(x)⁢log⁢p˜(x)q⁡(x;θ)⁢d⁢x(2)=∫p˜(x)⁢log⁢p˜(x)⁢d⁢x-p˜(x)⁢log⁢q⁡(x;θ)⁢d⁢x(3)=C-∫p˜(x)⁢q⁡(x;θ)q⁡(x;θ)⁢log⁢q⁡(x;θ)⁢d⁢x(4)=C-𝔼q⁡(x)[p˜(x)q⁡(x)⁢log⁢q⁡(x;θ)](5)≈C-∑i⁢w¯(i)(x(i);θ),(6)
where x(i)˜q(x; θ) and

w¯(i)=w(i)∑i⁢w(i);w(i)=p˜(x)q⁡(x;θ),
and C is a constant not depending on θ.

Applicant notes that the first term in equation (6) does not depend on the parameters θ; hence minimizing the KL divergence is equivalent to maximizing Σiw(i)log q(x(i); θ). The approaches have indicated that forward KL divergence may perform well in image rendering application.

Applicant submits that forward KL divergence does not promote sufficiently efficient sampling. For example, the KL divergence may be modelled by the Renyi divergence between the distributions, defined as:

Dα(q⁡(x)⁢p⁡(x))=1α-1⁢log⁢∫p⁡(x)α⁢q⁡(x)1-α⁢d⁢x(7)as α→1.

As is shown in eq. (7), the KL forward divergence does not provide a sufficiently severe penalty for any x where density of the target distribution p(x), is greater than the density of the importance distribution q(x), which discourages efficient sampling where samples have said characteristics.

Applicant proposes an approach that encourages the proposal distribution to explore the support of the target distribution as much as possible. Two alternate approaches are described; a first approach that determines error based on the Chi-Square distance between the target distribution and the importance distribution (which is equivalent to minimizing the variance of importance weights), and a second approach which further defines the error based on the KL divergence and selectively weighting said KL divergence.

By configuring error estimator112to determine loss as a function of the Chi-Square distance between the target distribution and the importance distribution, the error estimator112is more likely to explore the support of p than forward KL divergence. For example, optimizing the aforementioned Renyi divergence of α=2, is equivalent to optimize the log of variance of importance weights (or Chi-Square distance between the target and importance distribution), as shown below:

X2(f;θ)=log⁢𝕍q[p˜(x)q⁡(x)].(8)=log⁡(𝔼q[(p˜(x)q⁡(x))2]-[𝔼q(p˜(x)q⁡(x))]2)(9)=log(𝔼q[(p˜(x)q⁡(x))2]⁢-(∫p~(x)⁢dx)2)︸Constant(10)=log⁢𝔼q[(p⁡(x)Z⁢q⁡(x))2]-Constant(11)=log⁢𝔼q[(p⁡(x)q⁡(x))2]-log⁡(Z2)-Constant(12)=log⁢𝔼q[(p⁡(x)q⁡(x))2]-Constant.(13)

Where the equivalence between Chi-Square distance and minimizing is shown as equation (10) can be interpreted as minimizing the variance of the normalizing constant's estimator (That is, one chooses ƒ to be

f⁡(x)=∫p˜(x)⁢d⁢x=Z,thenlog⁢𝕍q[f⁡(x)⁢p˜(x)q⁡(x)]=log⁢𝕍q[Z⁢p˜(x)q⁡(x)]=log⁢𝕍q[p˜(x)q⁡(x)]-C,where⁢C⁢is⁢a⁢constant⁢).

This amounts to finding proposal distribution that gives the best estimate of the normalizing constant.

Applicant submits that, for α=2, which is a subset of instances that include an α larger α>1, the error estimator112penalizes more severely for any x such that p(x)>q(x).

Applicant proposes that the error estimator112operate according to the following function:

O𝒳2(F;θ)=∑jρ⁢O𝒳2(fj;θ)(14)
where each

O𝒳2(fj;θ)=𝕍[fj(x)⁢w]≈1N⁢∑i⁢(fj(x(i))⁢p˜(x(i))q(x(i))2.
where hereinafter

(x(i))q(x(i)
will be referred to as the importance weight, or importance value and denoted by w.

For each ƒj(where j is an integer) is defined as follows, ƒ0(x(i))=1 and

f1(x(i))=1D⁢∑d=1D⁢xd(i).

For

j>1⁢fj(x(i))=1D⁢∑d=1D⁢(xd(i)-x¯d)j
where the sample x(i)=(x1(i), x2(i), . . . , xD(i))˜q andx=(x1, . . . ,xd) is the empirical sample mean of proposal distribution.

Applicant's proposed approach avoids the technical challenge of past approaches in that different weights are assigned to different importance weights. For example, one can think of equation (8) as a special case ofq[ƒ(x)w] when ƒ(x)=1. That is, for each of x˜q(x), the penalty assigned to individual importance weights is constant (independent of whether the sample comes from the tail or central regions of q; seeFIG.10C).

The proposed error function is an error lossq[ƒ(x)w] in which ƒ(·) assigns each importance weight w a unique penalty based on the location of the value relative to the importance distribution geometry and encourages the importance distribution to match target distribution on multiple moment estimations (normalizing constant, mean, variance, skewness and kurtosis). For example, the i-th importance weight is assigned a positive scalar penalty

(1D⁢∑d=1D⁢(xd(i)-x_d)j)2
in the loss function, and samples x(i)that are further away from the empirical mean of importance distribution (i.e., samples in the tail regions of proposal) will be penalized based on their distance from the empirical mean of target distribution in the polynomial order (seeFIG.10C) and the gradient of

p⁡(x(i))q⁡(x(i))
is scaled up. In this way, the scalar function ƒjallows the approach to match the j-th of centered moment of the importance distribution with the target distribution.

The proposed error function strongly encourages space exploration and overcomes the pitfalls from the tail underestimation by incorporating a penalty based on the samples location under the importance distribution. Moreover, since poor estimation often results from the infinite variance of importance weights attributed to the tail underestimation, matching the centered moments (rather than un-centered ones) is important to reduce the possibility of tails of q falling below tails of p.

In example embodiments, other moment estimators are used to define the degree of the scalar penalty assigned to the sample based on its location in the importance distribution.

Applicant further proposes another method of operating the error estimator112. While the first approach promoted more efficient sampling, it may not necessarily lead to desirable performance on higher order estimations. Estimating error based solely on Chi-Square may encounter issues where the dimensionality of the target distribution becomes higher, as the performance of an error estimator based on Chi-Square distance becomes similar to the forward KL divergence, in the sense of not being able to explore the support of target distribution.

Alternatively stated, the first proposed approach may, in example embodiments, encourage aggressive exploration of the support of importance distribution, but possibly inhibit fine-tuning of the importance distribution to the local regions of target distribution.

Applicants propose to balance the trade-off between the competing goals of aggressive exploration and fine-tuning through the use of a linear annealing schedule between equation (14) and by incorporating KL divergence (more localized method as explained above) with the following formula:
O(β,θ)=βOx2(F;θ)+(1−β)OKL(θ)  (15)

The annealing schedule (i.e., β and (1−β) may be responsive to the different operating states of the MC estimator104, where in the first state, the MC estimator starts to update the model parameters of the sample generator106through training of the normalizing flow by optimizing Chi-Square distance (i.e., β is large or close to one) to take advantage of the greater stability to optimize and, in a second state, gradually shift the approach to error estimation towards the KL loss (i.e., β is small or close to zero) once the error estimator112becomes more stable in order to capture the tail behaviour.

Samples which are used in the MC Estimator104are generated by sample generator106.

Sample generator106includes an importance distribution model108A and a sample transformer108B. The Sample generator106is a computer representation of an importance distribution and a mapping function, and is alternatively referred to as the proposal distribution model data object.

The sample transformer108B is configured to draw samples from the importance distribution model108A and transform them into samples of the target distribution (or estimates of what would be target distribution samples) based on the trainable plurality of parameters.

In example embodiments, the sample generator106, which is representative of normalizing flow architecture, is configured to include one or more parameters to transform a random variable, z∈Dfrom the importance distribution model108A, through bijective mappings ƒθparameterized by θ (i.e., the sample transformer108B) into samples from the target distribution.

The sample transformer108B, given the probability density function of the variable z, represented as pz(z) (e.g., the function defining the base distribution of the importance distribution model108A), generates the probability of the transformed variable x=ƒθ(z), given by the multivariate change-of-variable formula:
qx(x)=qz(ƒθ−1(x))|degJ(x)|  (16)
where J is the Jacobian matrix of ƒθ−1.

Thus the density of x (a sample observation) is determined through a combination of the density of the importance distribution pzcomposed with the mapping ƒθand the local change of volume induced by each transformation in the form of the absolute value of the determinant of the Jacobian.

Hereinafter, for ease of reference, q(·; θ) will represent the sample generator106with the parameters θ. Applicants denote this model as FAIS (flow adaptive importance sampling).

As will be and has been described herein, the parameters θ of importance distribution model108A and sample transformer108B are iteratively updated via a gradient-based method and further based on a MC estimate evaluated with samples generated by sample generator106. Updating the parameters improves the performance of the importance distribution model108A and the sample transformer108B while computing MC estimates simultaneously.

In some embodiments, the proposed approach mitigates the technical challenge of using having a light tailed distribution as a base distribution, which has difficulty modelling heavy tailed distributions, by using a Student-t as the base distribution in the sample generator106. Applicants, during experimentation, observed that it is formidable challenge for a sample generator106with light-tailed base distribution (such as a normal distribution) to learn heavy-tailed target distributions. The sample generator106, with a Student-t as the base distribution, more accurately models heavy tailed distributions by learning the degree of freedom of the Student-t distribution based on the training data. In example embodiments, one of a variety of heavily tailed distributions may be used in place of a student-t distribution, such as Half-t distribution, Half-Cauchy and log-Normal distributions.

For example, when the degree of freedom of the Student-t base distribution approaches 1, the distribution converges to the Cauchy distribution with heavy tails. As the degree of freedom goes to infinity, the distribution approaches a Gaussian distribution. This allows the sample generator106, with a Student-t as the base distribution, to effectively capture a wide family of heavy-tailed distributions.

Returning again toFIG.1B, the method1000of MC estimation including updating the plurality of parameters of the machine learning architecture110includes:

At block1002A, the method includes randomly generating samples from an importance distribution defined by a plurality of parameters.

At block1002B, the method includes mapping, using a bijective mapping function defined by the plurality of parameters, the generated samples to a target distribution. The result of the sampling includes a data object which can represent an expected sample of the target distribution. For example, a sample could be a data structure or object storing a value representative of an interest rate value, a rate of infection, etc.

At block1006, the method includes determining a plurality of moment penalty terms (alternatively referred to as scalar penalty terms). For example, the plurality of moment penalty terms may be the moment penalty terms for a batch of samples, as described herein, and dependent upon the estimated mean of the samples of the batch.

At block1008, the method includes determining a loss with the loss function (e.g., equations 14 or 15). The determined loss is subsequently used to update the plurality of parameters.

To summarize, in a first state of operation, the machine learning architecture100is configured to generating samples at block1002, based on blocks1002A and1002B, determine importance weights based on the generated samples at block1004and determine moment penalty terms based on block1006. The machine learning architecture100is further configured to determine a loss based on the importance weights and the moment penalty terms, and update the plurality of parameters that define the bijective mapping function and the importance distribution based on said loss. In this first state, the importance weights are discarded.

In a second state of operation, the importance weights are further utilized to generate the MC estimates at block1010. For example, a MC estimates could be a data value in an output data structure, which data value could be representative of an expected price of an instrument, an expected amount of individuals to be infected, and so forth.

An algorithm representation of method1000is shown below:

Algorithm 1:Input : T := total number of training steps; O(Θ) := loss function; {tilde over (p)}(x) :=target distribution; q(x; Θ0) := normalizing flow parameterized byΘ0= {θ0; df0} where df0is the initial degree of freedom and θ is theflow parameters; M ∈ N number of exploratory stepsS ← { }; W ← { };for t := 0, . . . T do|Draw⁢{xt(i)}i=1S~q⁡(x;θt)⁢and⁢compute⁢weights⁢wt(i)←p~(x)q⁡(x;Θt)⁢∀i|Compute⁢loss⁢function⁢O(Θ)⁢and⁢gradient⁢∂O∂Θ| Update Θt+1using Adam| if t > M then| | Save the samples and weights S ← S ∪ {xt(i)}; W ← W ∪ {wt(i)}| endendreturn S, W

Referring now toFIG.1C, a schematic diagram of parallel sample generation with the system ofFIG.1A, according to example embodiments is shown.

Processor102, shown as a series of processors102A,102B, and102C (or a series of processing threads) which access or retrieve the importance distribution modeller108A and the sample transformer108B of sample generator106.

Each of processors102A,102B, and102C generate a sample (respectively sample observations154A,154B, and154C) independent of one another based on the current plurality of parameters defining the importance distribution modeller108A. Each of processors102A,102B, and102C generate a weight value for the respective sample (e.g., weight values156A,156B, and156C for sample observations154A,154B, and154C). In example embodiments, some processors are configured to continually generate samples, and other processors continually generate weights based on retrieved samples, all operating in parallel. Each processor may further be configured to generate, via error estimator112, a gradient value associated with the respective sample observation, wherein the aggregate gradient value of the batch of samples (i.e., the samples generated prior to updating the plurality of parameters) is used to update the plurality of parameters. For clarity, the gradient value associated with the respective sample observation is the evaluation of the gradient, for example, of equation 14 for the particular sample observation.

As presented herein, Applicants empirically demonstrate that the proposed approaches not only improves the range of distributions that can be accurately modelled with importance sampling, but also leads to improvements in at MC estimation and maximum likelihood.

In example embodiments, the operating characteristics of the proposed approach are measured with one or more diagnostics to determine whether the importance distribution matches the target distribution. An implicit advantage of implementing the sample generator106as a normalized flow architecture is that said architecture admits some natural diagnostics for assessing whether the proposal is matched up with target distribution.

For example, matching can be done by drawing some test samples x˜p(x), and then performing inverse function z=ƒ−1(x). If zqz, there is evidence that the proposal has fitted target poorly.

In the context of importance sampling, one cannot generate samples from target distribution. However, the system can use a sub-sampling method to assess the learning performance without accessing the samples of {tilde over (p)}.

FIGS.2A-2Care scatter plots that show how inversely transformed sub-samples behave based on how well q approximates the target distribution, shown as an example inFIG.2A.

The model inFIG.2Bapproximates the target distribution well and its inversely transformed sub-samples are similar to the base distribution. On the other hand, the model inFIG.2Cpoorly fits the target, and its inversely transformed sub-samples don't match Student-t distribution for any degrees.

Experiments

In example embodiments, the operating characteristics of the proposed approach are measured by measuring the error estimated by error estimator112against the total number of calls to {tilde over (p)}(x). The number of calls to this function may translate into the number of samples drawn from the importance distribution modeller108A. Measuring according to a number of samples drawn may be a good indicator of performance as, since all the MC methods are designed to be the asymptotic convergent, Applicants only interest is in the empirical convergent speed of the estimation. Measuring performance by the total number of samples required to achieve a certain error level places emphasis of the computing resources required to achieve a particular rate of error, as compared to the resource independent metric of error.

Assessing the error level in reference to the computational time, as compared to computational resources, may not be relevant, since for instance, Hamiltonian MC (HMC) is essentially a sequential method and inherently requires more time compared to annealed importance sampling (AIS) can be in parallel.

Moreover, comparison using the number of parameters of the model is not an appropriate choice, since HMC has only one parameter (step-size) and Stein Variational Gradient Descent (SVGD) are non-parametric. Moreover, Applicants do not use the effective sample size (ESS) as performance measure since HMC and (importance sampling) IS interpret ESS quite differently. Moreover, it is possible to have large ESS for models that are unable to estimate the second moment accurately.

The following section describes particulars of the Applicant's experimental testing:

RealNVP (Dinh et al., 2017) was used as the importance distribution modeller108A, and the sample transformer108B. The proposed method may be required to efficiently generate samples and perform density estimation, and RealNVP allows Applicants to do both sample estimation and density estimation in a single pass in parallel.

Many of Applicants' approaches to which the proposed approaches are compared are non-parametric models such as SVGD, AIS or methods such as HMC which have only one parameter (step size) while the sample generator106can contain hundreds of parameters. As a result, Applicants ignore the model complexity as performance measure (number of the parameters).

Methods such as HMC, which require evaluation of the exact gradient of target distribution {tilde over (p)}(x) for each sample, are treated in the experiments as requiring one additional call of target distribution because the additional evaluation provides more information about each sample.

Applicants have chosen six types of target distributions that are popularly used in the literature to study efficiency of sampling methods. These distributions include Gaussian mixture models, and some embodiments on estimating the marginal likelihood of Restricted Boltzmann Machine and variational Auto-Encoder. Both cases are zero order moment estimation and are generally easier than the other problems.

In the experimental testing, Applicants' method tackled various technical challenges described herein and outperformed the comparison approaches. The experimental results were robust to the hyper-parameters of MC estimator106(that is, Applicants use the same set of hyper-parameters across all problems without fine-tuning; the only exception made here is that a different network structure for the mixture model).

The Applicants experiments indicated that, even in high dimension, a small batch size (S=100≤D) at each iteration is sufficient to generate importance distributions that allow for acurate modelling of target distributions for many problems (with only the mixture model as exception). Having a small training batch size is crucial for deep learning models, because it allows one to build deeper network without too much concerns of computational complexity.

For first, second, third and fourth moments estimation, Applicants use the log Mean squared error (Log-MSE) for error measure. Given the true value and its estimator ƒ, {circumflex over (ƒ)}∈D, the error is defined as:

log⁢(1D⁢∑d⁢(f^d-fd)2).
For normalizing constant estimation, the error is defined as: log(|log(Z)−log({circumflex over (Z)})|).

The interpolation starts with β=1 and linearly decreasies down to some values close to 0 as the training process goes on. But β never goes to zero so that the exploration is kept as part of training. The polynomial order ρ in (equation (13)) was set to 4 during experimentation.

In this disclosure, Applicants show the results of 30 independent experiments, and report the averaged error in said experiments. Note that by doing this, the high standard deviation of estimator implicitly can be penalized, as discussed in respect of standard deviation below.

For the experimentation, the Applicants divided training into two stages which required the MC estimator104to operate in two states. In the first stage, the MC estimator104operates in a first state, reflective of exploration exploration. During the first stage, the importance distribution modeller108A can be poorly fit to the target distribution, and most of the samples generated by sample generator106can lead to arbitrarily bad MC estimations. During the first stage, MC estimator104is configured to discard samples generated by sample generator106when computing the MC estimation, but said samples are used for updating model parameters with error estimator112(which is similar to the burning-in stage of HMC/MCMC). During a second stage, the stage of exploitation, the MC estimator104operates according to a second state, wherein the samples are used both for updating model parameters with error estimator112, and for generating MC estimates.

In the experiments, the same two stage approach was used for other importance sampling approaches.

In these experiments, Applicants take the samples generated from the most recent 5×104or 105number of target evaluations as MC estimation samples.

For all error curves in the experiments, Applicants proposed approach in example embodiments exhibits a common phenomena: as more and more exploit samples are generated for MC estimation, the error curve of the sample generator106is relatively smoothly decreasing much faster than other methods. This may show that the proposed exploratory stage has successively approximated the target distribution quite well. The exploit samples now are safe to be used for the construction of reliable estimators.

10D Mixture Experiment

FIGS.3A-3Crelate to experiments where the target distribution was a 10 dimensional Bi-modes diagonal Gaussian distribution with variance one defined by

p~(x)=12⁢N⁡(x⁢❘"\[LeftBracketingBar]"μ1,I)+12⁢N⁡(x⁢❘"\[LeftBracketingBar]"μ2,I)⁢where⁢μ1=(2,…,2),μ2=(-2,…,-2).

FIG.3Ashows a visual comparison of the target distribution in a top row, and the learned importance distribution, in the bottom row, for three separate target distributions. The importance distribution was trained using 2D mixture models to model the complicated geometry.

FIG.3Bshows a histogram plot with the importance weights corresponding to samples at each iterations in log-space.FIG.3Bshows that as the training goes on, the variance of weights goes small and converges to zero. This indicates that Applicants importance distribution matches up with target distribution.

FIG.3Cshows an exaggerated visualization to visually represent the Applicants approach where, for the drawn samples x such that the target distribution304, represented by p(x). disagrees with the importance distribution302, represented as q(x), where x from the tail region of q(x) (e.g., samples308and306), will have a larger penalty assigned to them compared to samples closer to a mean of a distribution (e.g., sample310). This helps Applicants q to “lift-up” its tails.

Table 1, shown below, describes the mixture distribution result and the numbers in the table are the log error of MC estimation based on the learned importance distribution.

(r)1-2 Algorithmp[x]p[x]log(Z)FAIS−8.46−6.62−2.90HMC1.39722.777NAAIS−6.102−4.80−1.93GAIS1.3852.772−0.350LAIS1.3922.83620.390SVGD−6.17−3.07NAstepsize = 0.5

FIGS.4A-4Bdescribe the measured MC error on first (FIG.4A) and second moments (FIG.4B) in these experiments. These two estimations inform of the ability of the proposed approach to model on multi-modes in comparison to comparison methods.

Low estimation errors from SVGD, AIS and FAIS shows that they are capable of capturing the multi-modes (the improvement of FAIS indicates that FAIS still generates more efficient samples beyond just capture modes) while the biased estimations from AHMC, GAPIS, LAIS show that samples from these methods are stuck in local modes.

FIG.4Cdemonstrates the error of the normalizing constant for each algorithm, with lines402,404,406, and408indicating the results for LAIS, AIS, FAIS, and GAPIS respectively.

15D Correlated Schwefel Distribution Experiment

Schwefel Distribution and its similar version RoughWell contains an infinite number of local minimums and local maximums that can easily fail many adaptive approaches. Applicants make the problem more challenging by injecting high correlations between each dimension with randomly initialized means. Define the target distribution to be:

p~(x)∝exp⁡(-12⁢(x-μ)⁢∑-1(x-μ)T-∑i⁢sin⁡(π⁢❘"\[LeftBracketingBar]"xi❘"\[RightBracketingBar]"0.01))
where μ is randomly initialized uniformly in range [−5,5] and Σ is the covariance matrix value 10 along the diagonal and 9.95 elsewhere.

Table 2, shown below, describes the Schwefel Distribution result and the number in the table is the log error of MC estimation:

(r)1-2Algorithmp[x]ovp[x]urtosisp[x]FAIS−3.45−1.25−3.50AHMC−1.634.600.11AIS−0.784.5920.40GAIS0.554.6030.05LAIS−0.484.598.54SVGD1.184.481.69stepsize = 0.1

FIGS.5A-5Cdescribe the measured mean (FIG.5A), co-variance matrix (by averaging the error of each entry from co-variance matrix) (FIG.5B) and fourth moments estimation (FIG.5C), with lines502,504,506,508,510, and512indicating the results for LAIS, AIS, FAIS, GAPIS, AMHC and SVGD, respectivelty. AIS (with HMC transition) and HMC fails on such task.

The “Rough” surface of the distribution increases the challenge in the simulation of Hamiltonian dynamic with discretization leapfrog; hence, the properties of Schwefel Distributions lead to a very low acceptance rate. FAIS, as its optimization is independent of the exact gradient of target distribution, is able to circumvent the ill-behaviours of sinusoids and explore the correlated space.

100D Banana Distribution Experiment

The proposed method of some embodiments is discussed in references to non-linear correlated distribution, or Banana Distributions (BD), in this section. The challenge of this distribution arises from its strongly twisted shape that leads to the skewness, and high variance. The first moment also presents challenges to estimate because of the asymmetric two tail regions and requires sufficient exploration on both sides tails. The target is defined as

p~(x)=N⁡((x1,x2-b⁡(x22-σ2),…,xD)⁢❘"\[LeftBracketingBar]"0,∑)
where b=1; the covariance matrix Σ has ones along the diagonal and 0.95 for the rest of the entries. The major challenge of this target distribution is to estimate the mean, variance and skewness of second coordinate and the variance of all coordinates.

FIGS.6A-6Eshow the measured variance error averaged over dimensions (FIG.6A), and also, the statistics of the second coordinate (first (FIG.6C), second (FIG.6D), third moments (FIG.6E)) because these four are the most challenging tasks as discussed above, with lines602,604,606,608,610, and612indicating the results for LAIS, AIS, FAIS, GAPIS, AHMC, and SVGD, respectivelty. The low error on first, second and third moment estimation indicates that FAIS is capable of modeling the two asymmetric tails of target distribution.

On the other hand, the high error of second and third moments estimation from baseline methods suggests the following possible phenomena: (1) the sampler is stuck on one sided-tail, (2) the sampler can only model the central region, but is unable to reach out to the tail region.FIG.6Bdemonstrates the error of the normalizing constant for each algorithm.

Table 3, shown below, shows the log error of MC estimation based on a Banana Distribution:

(r)1-2Algorithmp[x2]p[x2]p[x]kewp[x2]log(Z)FAIS−10.66−6.842−9.56−7.958−1.778AHMC−1.2680.965−0.5511.714NAAIS−0.971.6361−0.9751.0850.092GAIS−0.651.911−0.6631.273−0.6381LAIS1.432.17310.0336.5826.648SVGD−1.810.901−1.584−1.570NAstepsize = 0.2

100D Funnel Distribution Experiment

Funnel Distribution (FD) is a special case of a hierarchical Bayesian model. It contains almost all pathological geometry such as varying curvatures, extreme outliers, asymmetry, etc. that can fail all types of sampling methods. Even worse, the high curvature region has high density with small volume such that many optimization methods either completely stuck in the low volume region and are unable to escape from it or traversing in the high volume region and become incapable of maneuvering into the low volume area.

Both cases lead to biased estimation. Define target distribution as:

p~(x)=∏d=2D⁢N⁡(xd⁢❘"\[LeftBracketingBar]"0,ex1)⁢N⁡(x1⁢❘"\[LeftBracketingBar]"0,9)

Note that the major challenge in Funnel Distribution lies in the estimation of the first order moment of first coordinates x1and the second order moment for all coordinates.

FIGS.7A-7D, show the measured first (FIG.7C) and second (FIG.7D) moments of first coordinate and the variance error averaged over dimension (FIG.7A) with lines702,704,706,708,710, and712indicating the results for LAIS, AIS, FAIS, GAPIS, AHMC, and SVGD respectivelty.

The high error of the estimation on first coordinate may be attributed to two possible scenarios: (1) the failure of collecting samples from funnel region (such as HMC, SVGD); (2) the samples collected are mostly from funnel region (trapped in the high curvature region; such as GAPIS method).

On the other hand, the averaged high error over dimensions on variance estimation attributed to the inefficiency of learning the tail region of target distribution (failed to traverse the low curvature region). Applicant notes here that FAIS is able to capture both low and high curvature regions relatively well compared against other baselines.FIG.7Bdemonstrates the error of the normalizing constant for each algorithm.

100D VAE Experiment

Quantitative analysis on generative model has been studied where AIS has been used for evaluating marginal likelihood of trained generative model. Applicants train a VAE model proposed in Appendix C.1 from on the MNIST dataset. The target distribution is:
{tilde over (p)}(z|x)∝p(x|z)p(z)
where

p⁡(z)⁢is⁢⁢N⁡(0,I)⁢and⁢p⁡(x⁢❘"\[LeftBracketingBar]"z)=exp⁡(∑1D⁢xi⁢log⁡(yi)+(1-xi)⁢log⁡(1-yi))
where x is the output from a two-layer decoder with ReLu activation function and y is the original image. In this case, marginal likelihood of x is the normalizing constant of above target distribution.

Table 4 shows the result, as the number in the table, being the log error of MC estimation based on a Funnel Distribution:

(r)1-2Algorithmp[x1]p[x1]p[x]log(Z)FAIS−4.07−0.343.37−3.5AHMC0.844.078.72NAAIS−2.153.868.94−0.33GAIS5.464.348.981.87LAIS−1.94.388.972.53SVGD1.13.948.98NAstepsize = 0.02
Restricted Boltzmann Machine

The estimate utilizing the partition function of RBM is discussed in this section. Let x∈Dbe the variable of interest and h∈{1, −1}10be hidden binary variables. Applicants randomly generated matrix W, and vector b, c from standard normal distribution. The target distribution is:

p~(x)∝∑h⁢exp⁡(xT⁢Wh+bT⁢x+cT⁢h-12⁢x2)

FIGS.8A and8Bshow the graphical results based on the above.

The performance on the marginal likelihood estimation for VAE is shown inFIG.8Aand for RBM is shown inFIG.8B, with lines802,804,806and808indicating the results for LAIS, AIS, FAIS, GAPIS respectivelty. FAIS can reduce the error below zero by utilizing the first 104number of {tilde over (p)}(x) evaluations while AIS requires 5×105.

APPENDIX

Loss Function Comparison

FIGS.9A-9Dassess the learning ability of flow model by measuring the log-likelihood, which describes the fit of q(x) to p(x). Lines902,904,906, and908ofFIGS.9A and9Bindicate the log-likelihood results for IAF with Guassian, RealNVP with Student-t, RealNVP with Guassian, and FFJORD with Guassian models, respectivelty.FIG.9Ashows 2D Pareto Banana Distribution MLE,FIG.9Bshows 2D Gaussian Banana Distribution MLE.

Lines920,922,924,926,928, and930ofFIG.9CandFIG.9Dshow the experimental log-likelihood of, respectively, approaches including: a Gaussian base distribution and a KL loss function, a Student's t base distribution and KL loss function, a Gaussian base distribution and a Chi Square loss function, a Student's t base distribution and a Chi Square loss function, a Gaussian base distribution and a polynomial loss function, and a Student's t base distribution and a polynomial loss function.

FIG.9Cshows 100D dimensional Funnel Distribution importance sampling.FIG.9Dshows 2D Gaussian Banana Distribution importance sampling.

FIG.9Ashows that IAF and FFJORD models are incapable of modelling heavy tailed distributions even in the 2 dimensional cases. RealNVP with Gaussian base performs better than IAF and FFJORD due to its explicitly exponential scaling; however, the occasional spikes indicates that it still fits the target poorly and such phenomena gets worse as the dimensionality increases. With Student-t based RealNVP, there may be a very smooth NLL curve.

FIG.9Bshows that when the target distribution is light tailed, Student-t based RealNVP can still perform equvalently good after few iterations. The earlier high loss is due to the fact that degree of freedom was initialized to be small at the beginning but the approach quickly adapts parameter to model the light tail.

FIG.9Cshows that a chart where, instead of asessing maximum likelihood on the generated data from target distribution, the importance sampling algorithm was executed and the likelihood was performed at each iteration. As in higher dimension, even with only one dimension is heavy tailed, Gaussian based RealNVP (e.g., lines920,922, and924) catastrophic failed to fit the target distribution while Student-t based distribution (lines926,928, and930) can fit well regardless of the choice of loss function.

FIG.9Dshows that, on light tailed distribution with importance sampling algorithm, that there isn't too much significant differences across all methods. Applicants Student-t based RealNVP can still successfully model the light tailed in this case.

FIGS.10A-10Cshow a comparison of Applicants optimization against KL, Chi-Square loss function. Error for each method (FIG.10Ashows the second moment estimation error on Funnel for lower dimensions;FIG.10Bshows the second moment estimation error on 100D Banana Distribution; andFIG.10Cshows the second moment estimation error on 100D Funnel Distribution) is averaged over 30 independent experiments. The estimation errors on low dimensional target distributions in which KL perform worse than Chi-Square; however, as dimension goes higher the difference is less significant while Applicants method is capable of outperforming on both low and high dimension problems (as shown inFIG.10BandFIG.10C).

Experiment Details

Flow Model

Applicants use the adam optimizer with learning rate 1e-3. The initial proposal distribution (base distribution) is Student-t distribution center at zero with diagonal variance of value 1. Applicants set the initial degree of freedom to be 64. 20 bijectors were used, each bijector consisting of 4 layers and each layer has 64 hidden units. The batch sample sized for each iteration is 100. Above three settings are for all experiments except for mixture model.

For the mixture model, Applicants use 5 bijectors, each bijector has 4 layers with 32 hidden units and the batch size is 10000 for each iterations.

AHMC

Given a fixed number of target evaluations, the first half of evaluations is for step-size adaption. The second half of evaluations is for HMC chain simulation. The initial samples is generated from N(0, I).

AIS

For the mixture model, Applicants set the batch size to be 1000 as this gives AIS the best estimation on Mixture. For the rest of the models, the batch size is fixed to 100, which is same as the proposed flow model approach. For all problems, each transition is done by simulating two iterations of HMC chain in which each iteration performs 2 number of leapfrog.

The initial proposal distribution is N(0, I).

SVGD

For each problem, Applicants ran SVGD for each choice of step size from {0.01,0.02,0.1,0.2,0.5}, and reported the best result for each model of them. Applicants used 1000 batch size for the mixture model. The batch size is set to be 100 (same as flow model) for the rest of the problems. The initial particles are generated from N(0, I).

GAPIS

Applicants used N(0,102) as the initial proposal distribution. Each batch size is 100 for all problems except mixture model (for mixture model, 10000 batch size is used). The gradient descent rate is 1−4.

LAIS

For each training iteration, Applicants performed 50 number of steps for MCMC warm-up, then with additional 200 steps of MCMC simulation. In the end Applicants collected 100 more samples from MCMC and used these samples to compute the empirical mean and co-variance matrix to update the proposal Gaussian distribution. The batch size setting is same as GAPIS and FAIS, that is 10000 samples for mixture model and 100 samples for other problems.

Standard Deviation

In this section, Applicants report another estimation accuracy measure. Consider the estimator {circumflex over (ƒ)}(i)results from the i-th experiment. The averaged estimator across all experiments defined as

f^_=1M⁢∑j⁢f^(j)
where M is the total number of experiments.{circumflex over (ƒ)}is a vector inD(in case that {circumflex over (ƒ)} is the estimator of normalizing constant, it's a constant). To analyze its accuracy, Applicants also look at the standard deviation over all experiments; which would be in the form of ({circumflex over (ƒ)}1±σ{circumflex over (ƒ)}1, . . . ,{circumflex over (ƒ)}D±σ{circumflex over (ƒ)}D)∈Dwhere{circumflex over (ƒ)}dand σ{circumflex over (ƒ)}dare the averaged estimated value and the standard deviation along d-th dimension respectively. It is possible to average over the dimensions; and define

μf^=1D⁢∑d=1D⁢f^_d⁢and⁢σf^⁢1D⁢∑d=1D⁢σf^d.
Table 5: Moisture Distribution Result: The Number in the Table is the Standard Deviation of the Estimator:

(r)1-2Algorithmp[x]p[x]log(Z)True050FAIS0.01 ± 0.015.00 ± 0.04−0.05 ± 0.07AHMC0.00 ± 2.20.99 ± 0.04NAAIS0.01 ± 0.085.00 ± 0.09−0.14 ± 0.05GAIS−1.99 ± 0.0051.00 ± 0.01−0.70 ± 0.03LAIS1.99 ± 0.230.89 ± 0.29−1.47 ± 0.05SVGD0.05 ± 0.074.78 ± 0.03NAstepsize = 0.5

Table 5: Note that although AHMC has seemingly achieved good results for the first moment estimation, it should be noted that the problem is equal weighted mixture model and the estimated value is averaged over many experiments in which each experiment stuck in one of the local modes; in this case, HMC's high standard deviation tells us that it does poorly on learning multi-modes.

Diagnostic

A natural (but not perfect) diagnostic for assessing the learning performance is to examine the disagreement between target distribution p and importance distribution q according to example embodiments is presented herein.

The method follows the re-sampling procedure with respect to the importance weights. By the end of training, a set of N samples is drawn, its importance weights are normalized to the proper probability distribution. Then, a set of M index is drawn with respect to weight distribution. The subset samples with respect to the index are retreived, and an inverse function is applied to invert this subset samples. Some multivariate statistics tests can be performed to assess the distribution of inverted samples whether the transformed samples follows the base distribution. If the inversely transformed sub-samples don't fit into base distribution, one can confidently claim q is poorly model p.

An example of the diagnostic method is shown below:

Algorithm 2: A simple diagnostic for flow based importance samplingmethodInput : A fully trained normalizing flow q(x; ΘT) parameterized by ΘT={θT; dfT}; N, M ∈ NDraw N number of samples {x(i)}i=1N~ q(x; θT) and compute weightsw(i)←p~(x)q⁡(x;ΘT)⁢∀iNormalize⁢the⁢weight⁢⁢w_(i)=w(i)∑i⁢w(i){x(i)}i=1M~ Multinomial(prob =w)Apply the inverse function z = fθT−1(x) on the re-samplesPerform Multivariate Student-T test with degree dfTto assess if {z}follow Student-T base distribution

TABLE 6Schwefel Distribution Result(r)1-2 Algorithmovp[x]urtosisp[x]TRUE9.95333FAIS9.65 ± 1.092.98 ± 0.36AHMC0.0001 ± 0.0012.59 ± 0.89AIS0.02 ± 0.069 · 106± 2 · 107GAIS0.002 ± 0.118 · 1013± 2 · 1014LAIS0.01 ± 0.041 · 103± 3 · 103SVGD24.97 ± 17.753.37 ± 2.79sstepsize = 0.1

Table 6: Note that AHMC's extremely small standard deviation and almost zero second moment estimation indicates, that adaptive method has adapted HMC's step-size so small in one of the local mode every run of experiment (another explanation could be that almost all samples get rejected).

TABLE 7Banana Distribution Result: the number in the table is the standarddeviation of the estimator:(r)1-2Algorithmp[x2]p[x2]p[x]kewp[x2]log(Z)TRUE−231.02−2.50FAIS−2.00 ± 0.012.99 ± 0.081.02 ± 0.01−2.54 ± 0.06−0.16 ± 0.01AHMC−2.10 ± 1.070.96 ± 1.070.27 ± 0.11−0.21 ± 0.22NAAIS−1.57 ± 0.500.69 ± 0.430.43 ± 0.09−0.77 ± 0.66−1.10 ± 0.02GAIS−1.27 ± 0.070.39 ± 0.080.33 ± 0.05−0.68 ± 0.16−0.53 ± 0.07LAIS−0.27 ± 1.620.03 ± 0.080.02 ± 0.07−1 · 104± 4 · 104−771 ± 48SVGD−1.60 ± 0.011.44 ± 0.040.58 ± 0.002−2.10 ± 0.13NAstepsize = 0.2

Table 7: Once again the HMC seemingly has good results on estimate the first moment of the second coordinate. However, large variance among the experiment runs indicates the instability of HMC.

TABLE 8Funnel Distribution Result(r)1-2Algorithmp[x1]p[x1]p[x]log(Z)TRUE098993FAIS0.11 ± 0.148.2 ± 0.5986.53 ± 5.7892.97 ± 0.04AHMC0.97 ± 2.951.24 ± 1.213.86 ± 24.16NAAIS−0.19 ± 0.332.09 ± 0.361.85 ± 0.8792.28 ± 0.06GAIS−15.36 ± 0.20.21 ± 0.070.002 ± 0.0086.48 ± 0.35LAIS0.42 ± 0.160.05 ± 0.050.62 ± 0.6180.34 ± 1.55SVGD−1.93 ± 1.121.69 ± 1.540.13 ± 0.16NAstepsize = 0.02

TABLE 9Restricted Boltzmann MachineAlgorithmlog(Z)TRUE220FAIS219.02 ± 0.04GAPIS160.07 ± 4.20LAIS174.04 ± 10.73AIS218.62 ± 0.07

TABLE 10Variational Auto-EncoderAlgorithmlog(Z)TRUE−76FAIS−77.42 ± 0.01GAPIS−84.01 ± 2.90LAIS−156 ± 10.77AIS−77.84 ± 0.03

InFIG.11, a system for computational estimation sampling from non-trivial probability distributions1100is shown, according to some embodiments.

An input x which is indicative of a target distribution is received at input receiver1102. In example embodiments, the target distribution is represented by a high dimensional integrals.

An importance distribution architecture1104is configured to receive the input to conduct importance sampling by drawing samples from an importance distribution, the importance sampling using a loss function adapted to minimize a variance of an importance sampling estimator for higher order moments.

In example embodiments, an importance distribution determination engine1106determines or selects the relevant importance distribution, an importance sampling estimator engine1108generates samples of a distribution based on the importance distribution, an importance sampling determination engine1110conducts importance sampling by reviewing samples from the importance sampling estimator engine1108, using a loss function adapted to minimize a variance of the importance sampling estimator1108samples for higher order moments.

The importance distribution architecture1104may further comprise a MC estimate generator1112configured to weigh the samples reviewed by the importance sampling determination engine1110by a ratio of a target posterior and the importance distribution to construct a MC estimate Y.

In example embodiments, the system1100may be configured to be a stand alone system. In some embodiments, system1100may be configured to connect to utilize the provided as a system interoperating upstream and downstream computing devices within a data center.

FIG.12is a schematic diagram of the computing device100, which may be a server or series of servers, which in some embodiments, may be a special purpose machine that is specifically configured for generating MC estimates, having, for example, specialized hardware components and/or software components, such as graphics cards or dedicated high performance processors adapted for parallel processing. The computing device100can be utilized as a computer server or computational mechanism for conducting and/or retrofitting MC estimates such that an improved distribution can be used for adaptive importance sampling, which can aid in establishing improved probability distributions for modelling distributions that have significant tails. This can aid the MC estimation from a computational perspective as it can improve accuracy for a given amount of available computing resources or processing capabilities.

As depicted, the computing device100includes at least one processor1202, memory1208, at least one I/O interface1206, and at least one network interface1208.

Processor1202may be an Intel or AMD x86 or x64, PowerPC, ARM processor, or the like. Memory1204may include a combination of computer memory that is located either internally or externally such as, for example, random-access memory (RAM), read-only memory (ROM), compact disc read-only memory (CDROM).

Each I/O interface1206enables computing device100to interconnect with one or more input devices, such as a keyboard, mouse, camera, touch screen and a microphone, or with one or more output devices such as a display screen and a speaker. A graphical interface can be provided which receives inputs and can be used to receive parameters and/or instruction sets to control retrofitting of computational models.

Each network interface1208enables computing device100to communicate with other components, to exchange data with other components, to access and connect to network resources, to serve applications, and perform other computing applications by connecting to a network (or multiple networks) capable of carrying data including the Internet, Ethernet, plain old telephone service (POTS) line, public switch telephone network (PSTN), integrated services digital network (ISDN), digital subscriber line (DSL), coaxial cable, fiber optics, satellite, mobile, wireless (e.g. Wi-Fi, WiMAX), SS7 signaling network, fixed line, local area network, wide area network, and others.

The network interface1208can connect to a message bus, for example, which may couple with downstream or upstream computational systems which may provide data sets for training and/or models for retrofit (e.g., data structures representing same).

Applicants, in example embodiments, propose a system where the importance distribution is adapted (according to the loss function) throughout the overall process of sampling.

Applicant notes that the described embodiments and examples are illustrative and non-limiting. Practical implementation of the features may incorporate a combination of some or all of the aspects, and features described herein should not be taken as indications of future or existing product plans. Applicant partakes in both foundational and applied research, and in some cases, the features described are developed on an exploratory basis.

The term “connected” or “coupled to” may include both direct coupling (in which two elements that are coupled to each other contact each other) and indirect coupling (in which at least one additional element is located between the two elements).

Although the embodiments have been described in detail, it should be understood that various changes, substitutions and alterations can be made herein without departing from the scope. Moreover, the scope of the present application is not intended to be limited to the particular embodiments of the process, machine, manufacture, composition of matter, means, methods and steps described in the specification.

As one of ordinary skill in the art will readily appreciate from the disclosure, processes, machines, manufacture, compositions of matter, means, methods, or steps, presently existing or later to be developed, that perform substantially the same function or achieve substantially the same result as the corresponding embodiments described herein may be utilized. Accordingly, the appended claims are intended to include within their scope such processes, machines, manufacture, compositions of matter, means, methods, or steps.

As can be understood, the examples described above and illustrated are intended to be exemplary only.