RareEventProbabilityEstimationContext Class

A RareEventProbabilityEstimator is executed by calling its Estimate method. This is a template method, in the sense that it defines the invariant parts of a Cross-Entropy program for rare event simulation. Such method relies on an instance of class RareEventProbabilityEstimationContext, which is passed as a parameter to it and specifies the primitive operations of the template method, i.e. those varying steps of the algorithm that depends on the problem under study.

Class RareEventProbabilityEstimationContext thoroughly defines the rare event of interest, and supplies the primitive operations as abstract methods, that its subclasses will override to provide the concrete behavior of the estimator.

A Cross-Entropy estimator is designed to evaluate the probability of rare events regarding the performance of complex systems. It is assumed that such probability must be evaluated with respect to a specific density function, member of a parametric family



where is a vector representing a possible state of the system, and is the set of allowable values for parameter .

Let be an extreme performance value, referred to as the threshold level, and consider the event defined as



or



where is the state space of the system and is the function returning the system performance at state . Let us assume that the probability of an event having form or must be evaluated under the specific parameter value , the so called nominal parameter: such probabilities can be evaluated by a Cross-Entropy estimator.

At instantiation, the constructor of a RareEventProbabilityEstimationContext object will receive information about the rare event under study by means of parameters representing , , and a constant stating if the event has the form or .

After construction, property ThresholdLevel represents , while can be inspected via property InitialParameter. Lastly, property RareEventPerformanceBoundedness signals that the event has the form if it evaluates to the constant Lower, or that it has the form if it evaluates to the constant Upper.

Implementing a context for rare event simulation

The Cross-Entropy method provides an iterative multi step procedure, where, at each iteration , the sampling step is executed in order to generate diverse candidate states of the system, sampled from a distribution characterized by the reference parameter of the iteration, say . Such sample is thus exploited in the updating step, in which a new level and a new reference parameter are identified to modify the distribution from which the samples will be obtained in the next iteration, in order to improve the probability of sampling relevant states, i.e. those states corresponding to the performance excesses of interest (See the documentation of class CrossEntropyProgram for a thorough discussion of the Cross-Entropy method).

When the Cross-Entropy method is applied in an rare event context, the estimating step is executed, in which the rare event probability is effectively estimated.

These steps must be implemented as follows.

Sampling step

Since class RareEventProbabilityEstimationContext derives from CrossEntropyContext, property StateDimension will return the dimension of the points in the Cross-Entropy samples. If a state of the system can be represented as a vector of length , as in , then should be returned. In addition, method PartialSample(Double, TupleInt32, Int32, RandomNumberGenerator, DoubleMatrix, Int32) must be overriden to determine how to draw the sample locally to a given thread when processing in parallel the sampling step.

Updating step

Class RareEventProbabilityEstimationContext overrides for you the methods required for ordering the system performances of the states sampled in the previous step, and for updating the iteration levels. However, to complete the implementation of the updating step, function must be defined via method Performance(DoubleMatrix), and method UpdateParameter(LinkedListDoubleMatrix, DoubleMatrix) also needs to be overridden. Given and , the method is expected to return the solution to the following program:



where



is the Likelihood ratio of to evaluated at , is the set of elite sample positions (row indexes of the matrix returned by method Sample, with being the -th row of such matrix and being its number of rows. Method UpdateParameter receives two parameters. The first is a LinkedListT of DoubleMatrix instances, representing the reference parameters applied in previous iterations. The instance returned by property Last corresponds to parameter . The second method parameter is a DoubleMatrix instance whose rows represent the elite points sampled in the current iteration.

Estimating step

The estimating step is executed after that the underlying Cross-Entropy program has converged. Given the optimal parameter , a Likelihood estimator based on the ratio is exploited by drawing a sample from density . Method GetLikelihoodRatio(DoubleMatrix, DoubleMatrix, DoubleMatrix) must be overridden in order to evaluate the ratio at a specific sample point for a given pair of nominal and reference parameters.

The following example is based on the Rare-Event Simulation Example proposed by Rubinstein and Kroese (Section 2.2.1)^[1] .

A Cross-Entropy estimator is exploited to analyze, under extreme conditions, the length of the shortest path from node to node in the undirected graph depicted in the following figure.

100%Shortest path from A to BRareEventExample.png

There are 5 edges in the graph. Their lengths are represented as independent and exponentially distributed random variables , with variable having a specific mean , for , forming the following nominal parameter:

In order to represent a possible state of such system, we hence need a vector of length 5:



and the performance of the system is the total length of the shortest path given , i.e.

The likelihood ratio of to the generic density is

At iteration , the parameter updating formula is, for ,



where .

In this example, we want to estimate the probability of observing a shortest path greater than or equal to a level set as . This is equivalent to define the target event as .

using System;
using System.Collections.Generic;
using Novacta.Analytics.Advanced;
using System.Linq;

namespace Novacta.Analytics.CodeExamples.Advanced
{
    public class CrossEntropyEstimatorExample0  
    {
        class RareShortestPathProbabilityEstimation
            : RareEventProbabilityEstimationContext
        {
            public RareShortestPathProbabilityEstimation()
                : base(
                    // There are 5 edges in the network under study.
                    // Hence we need a vector of length 5 to represent
                    // the state of the system, i.e. the observed lengths 
                    // corresponding to such edges.
                    stateDimension: 5,
                    // Set the threshold level and the rare event 
                    // boundedness in order to target the probability 
                    // of the event {2.0 <= Performance(x)}.
                    thresholdLevel: 2.0,
                    rareEventPerformanceBoundedness:
                        RareEventPerformanceBoundedness.Lower,
                    // Define the nominal parameter of interest.
                    initialParameter: DoubleMatrix.Dense(1, 5,
                        new double[] { 0.25, 0.4, 0.1, 0.3, 0.2 })
                      )
            {
            }

            // Define the performance function of the 
            // system under study.
            protected override double Performance(DoubleMatrix x)
            {
                // Compute the lengths of the possible 
                // paths when the state of the system is x.
                DoubleMatrix paths = DoubleMatrix.Dense(4, 1);
                paths[0] = x[0] + x[3];
                paths[1] = x[0] + x[2] + x[4];
                paths[2] = x[1] + x[4];
                paths[3] = x[1] + x[2] + x[3];

                // Compute the shortest path.
                var indexValuePair = Stat.Min(paths);

                // Return the shortest path.
                return indexValuePair.Value;
            }

            // Set how to sample system states
            // given the specified parameter. 
            protected override void PartialSample(
                double[] destinationArray,
                Tuple<int, int> sampleSubsetRange,
                RandomNumberGenerator randomNumberGenerator,
                DoubleMatrix parameter,
                int sampleSize)
            {
                // Must be Item1 included, Item2 excluded
                int subSampleSize = sampleSubsetRange.Item2 - sampleSubsetRange.Item1;
                int leadingDimension = Convert.ToInt32(sampleSize);
                for (int j = 0; j < this.StateDimension; j++)
                {
                    var distribution = new ExponentialDistribution(1.0 / parameter[j])
                    {
                        RandomNumberGenerator = randomNumberGenerator
                    };
                    distribution.Sample(
                         subSampleSize,
                         destinationArray, j * leadingDimension + sampleSubsetRange.Item1);
                }
            }

            // Define the Likelihood ratio for the 
            // problem under study.
            protected override double GetLikelihoodRatio(
                DoubleMatrix samplePoint,
                DoubleMatrix nominalParameter,
                DoubleMatrix referenceParameter)
            {
                var x = samplePoint;
                var u = nominalParameter;
                var v = referenceParameter;
                var prod = 1.0;
                var sum = 0.0;
                for (int j = 0; j < x.Count; j++)
                {
                    prod *= v[j] / u[j];
                    sum += x[j] * (1.0 / u[j] - 1.0 / v[j]);
                }

                return Math.Exp(-sum) * prod;
            }

            // Set how to update the parameter via 
            // the elite sample.
            protected override DoubleMatrix UpdateParameter(
                LinkedList<DoubleMatrix> parameters,
                DoubleMatrix eliteSample)
            {
                var nominalParameter = parameters.First.Value;
                var referenceParameter = parameters.Last.Value;

                var ratios = DoubleMatrix.Dense(1, eliteSample.NumberOfRows);

                var u = nominalParameter;
                var w = referenceParameter;
                for (int i = 0; i < eliteSample.NumberOfRows; i++)
                {
                    var x = eliteSample[i, ":"];
                    ratios[i] = this.GetLikelihoodRatio(x, u, w);
                }

                var newParameter = (ratios * eliteSample);
                var sumOfRatios = Stat.Sum(ratios);
                newParameter.InPlaceApply(d => d / sumOfRatios);

                return newParameter;
            }
        }

        public void Main()
        {
            // Create the context.
            var context = new RareShortestPathProbabilityEstimation();

            // Create the estimator.
            var estimator = new RareEventProbabilityEstimator()
            {
                PerformanceEvaluationParallelOptions = { MaxDegreeOfParallelism = 1 },
                SampleGenerationParallelOptions = { MaxDegreeOfParallelism = 1 }
            };

            // Set estimation parameters.
            double rarity = 0.1;
            int sampleSize = 1000;
            int finalSampleSize = 10000;

            // Solve the problem.
            var results = estimator.Estimate(
                context,
                rarity,
                sampleSize,
                finalSampleSize);

            // Show the results.
            Console.WriteLine("Under the nominal parameter:");
            Console.WriteLine(context.InitialParameter);
            Console.WriteLine("the estimated probability of observing");
            Console.WriteLine("a shortest path greater than 2.0 is:");
            Console.WriteLine(results.RareEventProbability);

            Console.WriteLine();
            Console.WriteLine("Details on iterations:");

            var info = DoubleMatrix.Dense(
                -1 + results.Parameters.Count,
                1 + results.Parameters.Last.Value.Count);

            info.SetColumnName(0, "Level");
            for (int j = 1; j < info.NumberOfColumns; j++)
            {
                info.SetColumnName(j, "Param" + (j - 1).ToString());
            }

            int i = 0;
            foreach (var level in results.Levels)
            {
                info[i++, 0] = level;
            }

            var referenceParameters = results.Parameters.Skip(1).ToList();
            var paramIndexes = IndexCollection.Range(1, info.NumberOfColumns - 1);
            for (i = 0; i < info.NumberOfRows; i++)
            {
                info[i, paramIndexes] = referenceParameters[i];
            }

            Console.WriteLine();
            Console.WriteLine(info);
        }
    }
}

// Executing method Main() produces the following output:
// 
// Under the nominal parameter:
// 0.25             0.4              0.1              0.3              0.2              
// 
// 
// the estimated probability of observing
// a shortest path greater than 2.0 is:
// 1.339307437308636E-05
// 
// Details on iterations:
// 
// [Level]          [Param0]         [Param1]         [Param2]         [Param3]         [Param4]         
// 0.5701086        0.483438676      0.736805982      0.116523615      0.480972597      0.351788383      
// 1.00101664       0.778124668      1.03296697       0.139379404      0.611880162      0.445196897      
// 1.43228893       1.12582143       1.36538794       0.135058159      0.663007916      0.532245242      
// 1.90265032       1.4492325        1.59888059       0.114423424      0.800780865      0.682784109      
// 2                1.87183229       1.86492053       0.115583952      0.913234886      0.802443085      
// 

Reference