The PythonRandomVector.getSample() method returns a 0-dimension sample

[mbaudin47]

What happened?

The QuantileOrderStatisticsBound class below implements a PythonRandomVector. More precisely, it implements the getRealization() method. The getSample() method is used in the script, but this produces an error: the returned sample has the correct size, but has a 0 dimension.

If we implement the getSample() method naively, everything works. Still, it seems to me that there should be a default implementation of getSample(), since that it just a loop over the required sample size.

How to reproduce the issue?

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):
        sample = []
        for i in range(resample_size):
            sample.append(self.getRealization())
        return sample
"""

# %%
x = ot.Uniform(0.0, 1.0)
t = 500 * x ** 5 + 1000
t.setDescription(["T"])
alpha = 0.95
xAlpha = t.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(
    t, alpha, sample_size, rank_index))
sample = random_vector.getSample(resample_size)
print(f"Mean = {sample.computeMean()[0] : .4f}")


produces:

IndexError: vector::_M_range_check: __n (which is 0) >= this->size() (which is 0)

Version

1.22

Operating System

Windows

Installation media

conda

Additional Context

No response

[jschueller]

there is a default implementation of getSample but it wants to use the value 0 you declared for the dimension:
super(QuantileOrderStatisticsBound, self).__init__(0)