We should be able to build a ThresholdEvent() on a RandomVector()

[mbaudin47]

What is the idea?

In the script below, we estimate the probability $p = \mathbb{P}\left(x_\alpha \leq X_{(59})\right)$ where $X$ is a random variable with uniform(-1,1) distribution, $x_\alpha$ is the quantile of $X$ at level $\alpha = 0.95$, { $X_1, ..., X_{59}$ } is a sample of $X$ with size 59 and $X_{(59)}$ is the largest observation in the sample. This is why we are interested by the event ${x_\alpha \leq X_{(59)}}$. To do this, we implement a Bernoulli random variable defined by:

$$ B = \begin{cases} 1 & \textrm{ if } x_\alpha \leq X_{(59)}\\ 0 & \textrm{ otherwise.} \end{cases} $$

It is straightforward to implement the associated PythonRandomVector. Based on this, we should be able to use the ProbabilitySimulationAlgorithm to estimate the probability based on a pre-defined tolerance on the accuracy of the probability estimate. Of course, we can still use the PythonRandomVector.getSample(1000).computeMean() method to estimate the probability, but we do not benefit from the features of the ProbabilitySimulationAlgorithm class:

• there is no iterative algorithm based on a pre-defined accuracy tolerance,
• there is no Result with confidence intervals, etc.

The problem is that the random vector is not of CompositeRandomVector type. Hence, we cannot create the ThresholdEvent required by the ProbabilitySimulationAlgorithm.

Why is this needed?

That feature is need when we want to check that some confidence interval of a level - $\alpha$ quantile has at least the required $\beta$ confidence level, using Monte-Carlo method.

Additional Context

The script:

# %%
import openturns as ot

# %%
class QuantileOrderStatisticsBound(ot.PythonRandomVector):
    def __init__(self, distribution, alpha, sample_size, rank_index):
        """
        Create a RandomVector of bound of a quantile from order statistics.

        Let x_alpha be the quantile of X of level alpha.
        Let X[k] be an order statistics of the
        random variable X. 

        Parameters
        ----------
        distribution : ot.Distribution
            The distribution of X.
        alpha : float, in [0, 1]
            The quantile level.
        sample_size : int
            The sample size.
        rank_index : int, in {1, ..., sample_size}
            The index of the observation in the sorted mathematical sample.
            The Python index is rank_index - 1 in {0, ..., sample_size - 1}.
        """
        super(QuantileOrderStatisticsBound, self).__init__(0)
        self.sample_size = sample_size
        if rank_index < 0:
            raise ValueError(f"The minimum rank {rank_index} is lower than 0.")
        if rank_index >= sample_size + 1:
            raise ValueError(
                f"The minimum rank {rank_index} is greater "
                f"than the sample size {self.sample_size}."
            )
        self.rank_index = rank_index
        self.distribution = distribution
        self.alpha = alpha
        self.xAlpha = distribution.computeQuantile(alpha)[0]

    def getRealization(self):
        """
        Generate a single observation of the event.

        Returns
        -------
        isOrderStatisticsBound : float, in {0, 1}
            Returns 1 if the event is observed.
            Returns 0 otherwise.
        """
        sample = self.distribution.getSample(self.sample_size)
        sortedSample = sample.sort()
        if (
            self.xAlpha <= sortedSample[self.rank_index - 1, 0]
        ):
            isOrderStatisticsBound = 1.0
        else:
            isOrderStatisticsBound = 0.0
        return [isOrderStatisticsBound]

    def getSample(self, resample_size):
        """
        Generate sample of the event.

        Parameters
        ----------
        resample_size: int
            The sample size
        
        Returns
        -------
        sample : ot.Sample(resample_size, 1)
            The sample of observations of the binary random variable.
        """
        sample = []
        for i in range(resample_size):
            sample.append(self.getRealization())
        return sample


# %%
distribution = ot.Uniform()
alpha = 0.95
xAlpha = distribution.computeQuantile(alpha)
print(f"Quantile at level {alpha:.4f} = {xAlpha[0]:.2f}")

# %%
sample_size = 59
resample_size = 10000
rank_index = 59

# %%
# Evalue la vraie probabilité
binomial = ot.Binomial(sample_size, alpha)
probability = binomial.computeCDF(rank_index - 1)
print(fr"P(x_\alpha <= X({rank_index})) = {probability:.4f}")

# %%
random_vector = ot.RandomVector(QuantileOrderStatisticsBound(
    distribution, alpha, sample_size, rank_index))
sample = random_vector.getSample(resample_size)
print(f"Mean = {sample.computeMean()[0] : .4f}")

# %%
event = ot.ThresholdEvent(random_vector, ot.Greater(), 0.5)

# %%
algo = ot.ProbabilitySimulationAlgorithm(event)
algo.setMaximumOuterSampling(150)
algo.setBlockSize(4)
algo.setMaximumCoefficientOfVariation(0.1)
# Perform the simulation
algo.run()
print('Probability estimate=%.6f' % algo.getResult().getProbabilityEstimate())

produces the error:

TypeError: InvalidArgumentException : Event can only be constructed from composite random vectors.

[jschueller]

similar to #2378