Documentation / Manuel utilisateur en Python :

As stated in Section II.2.2, these information are now stored in vectors because of the new possible attribute nature. In the case of constant-size vectors whose dimension is , the attribute-statistical vectors are filled with the statistical information considering every elements of the input vector independent from the others. This results in attribute-statistical vectors of size , the extra element being the statistical information computed over the complete vector. Here is an example of computeStatistic use with the tdstest.dat file already shown in Section II.3.4:

"""
Example of statistical estimation on vector attributes
"""
from rootlogon import DataServer

tdsop = DataServer.TDataServer("foo", "poet")
tdsop.fileDataRead("tdstest.dat")

# Considering every element of a vector independent from the others
tdsop.computeStatistic("x")
px = tdsop.getAttribute("x")

print("min(x[0])="+str(px.getMinimum(0))+";  max(x[0])="+str(px.getMaximum(0))
      + ";  mean(x[0])="+str(px.getMean(0))+";  std(x[0])="+str(px.getStd(0)))
print("min(x[1])="+str(px.getMinimum(1))+";  max(x[1])="+str(px.getMaximum(1))
      + ";  mean(x[1])="+str(px.getMean(1))+";  std(x[1])="+str(px.getStd(1)))
print("min(x[2])="+str(px.getMinimum(2))+";  max(x[2])="+str(px.getMaximum(2))
      + ";  mean(x[2])="+str(px.getMean(2))+";  std(x[2])="+str(px.getStd(2)))
print("min(xtot)="+str(px.getMinimum(3))+";  max(xtot)="+str(px.getMaximum(3))
      + ";  mean(xtot)="+str(px.getMean(3))+";  std(xtot)="+str(px.getStd(3)))

# Statistic for a single realisation of a vector, not considering other events
tdsop.addAttribute("Min_x", "Min$(x)")
tdsop.addAttribute("Max_x", "Max$(x)")
tdsop.addAttribute("Mean_x", "Sum$(x)/Length$(x)")

tdsop.scan("x:Min_x:Max_x:Mean_x", "", "colsize=5 col=2:::6")

The first computation is filling the vector of statistical elementary in the concerned attribute . The first, second and third cout line in the previous piece of code are dumping the statistical characteristics respectively for the first, second and third element of the vector. The fourth one is giving the main characteristics considering the complete vector and all the entries. The results of this example are shown below:

min(x[0])= 1;  max(x[0])= 7;  mean(x[0])= 3;  std(x[0])= 3.4641
min(x[1])= 2;  max(x[1])= 8;  mean(x[1])= 4.66667;  std(x[1])= 3.05505
min(x[2])= 3;  max(x[2])= 9;  mean(x[2])= 6.66667;  std(x[2])= 3.21455
min(xtot)= 1;  max(xtot)= 9;  mean(xtot)= 4.77778;  std(xtot)= 3.23179

This implementation has been chosen as ROOT offers access to another way of computing these notions if one wants to consider every element of a vector, assuming that every event is now independent from the others. Indeed it is possible to get the minimum, maximum and mean of a vector on an event-by-event basis by introducing a new attribute with a formula, as done in Section II.3.4. This is the second part of the code shown in the box above (using specific function from ROOT, that needs the sign "$" to be recognised). The results are shown below:

*****************************************************
*    Row   * Instance *  x * Min_x * Max_x * Mean_x *
*****************************************************
*        0 *        0 *  1 *     1 *     3 *      2 *
*        0 *        1 *  2 *     1 *     3 *      2 *
*        0 *        2 *  3 *     1 *     3 *      2 *
*        1 *        0 *  7 *     7 *     9 *      8 *
*        1 *        1 *  8 *     7 *     9 *      8 *
*        1 *        2 *  9 *     7 *     9 *      8 *
*        2 *        0 *  1 *     1 *     8 * 4.3333 *
*        2 *        1 *  4 *     1 *     8 * 4.3333 *
*        2 *        2 *  8 *     1 *     8 * 4.3333 *
*****************************************************

For a given probability , the corresponding quantile is given by:

where is the k-Th smallest value of the attribute set-of-value (whose size is ). The way is computed is discussed later on, as a parameter of the functions.

The implementation and principle has slightly changed in order to be able to cope with vectors (even though the previous logic has been kept for consistency and backward compatibility). Let's start with an example of the way it was done with the two main methods whose name are the same but differ by their signature.

from ctypes import c_double
import numpy as np

aproba = c_double(0.5)
aquant = c_double(0)
tdsGeyser.computeQuantile("x2", aproba, aquant) 

Proba = np.array([0.05,0.95], dtype=np.float64)
Quant = np.array([0.0,0.0], dtype=np.float64) 
tdsGeyser.computeQuantile("x2", 2, Proba, Quant)

print("Quant[0]=%f; Quant[1]=%f" % (Quant[0],Quant[1])) 	

This implementation has been slightly modified for two reasons: to adapt the method to the case of vectors and to store easily the results and prevent from recomputing already existing results. Even though the previous behaviour is still correct, the information is now stored in the attribute itself, as a vector of map. For every element of a vector, a map of format map<double,double> is created: the first double is the key, meaning the value of probability provided by the user, while the second double is the results. It is now highly recommended to use the method of the TAttribute, that gives access to these maps for two reasons: the results provided by the methods detailed previously are only correct for the last element of a vector, and the vector of map just discussed here is cleared as soon as the general selection is modified (as for the elementary statistical-vectors discussed in Section II.4.3). The next example uses the following input file, named aTDSWithVectors.dat:

#NAME: cho
#COLUMN_NAMES: x|rank
#COLUMN_TYPES: D|V

0 0,1
1 2,3
2 4,5
3 6,7
4 8,9

From this file, the following code (that can be find in Section XIV.3.14) shows the different methods created in the attribute class in order for the user to get back the computed values:

"""
Example of quantile estimation for many values at once
"""
from sys import stdout
from ctypes import c_int, c_double
import numpy as np
from rootlogon import ROOT, DataServer

tdsvec = DataServer.TDataServer("foo", "bar")
tdsvec.fileDataRead("aTDSWithVectors.dat")

probas = np.array([0.2, 0.6, 0.8], 'd')
quants = np.array(len(probas)*[0.0], 'd')
tdsvec.computeQuantile("rank", len(probas), probas, quants)

prank = tdsvec.getAttribute("rank")
nbquant = c_int(0)
prank.getQuantilesSize(nbquant)  # (1)
print("nbquant = " + str(nbquant.value))

aproba = c_double(0.8)
aquant = c_double(0)
prank.getQuantile(aproba, aquant)  # (2)
print("aproba = " + str(aproba.value) + ", aquant = " + str(aquant.value))

theproba = np.array(nbquant.value*[0.0], 'd')
thequant = np.array(nbquant.value*[0.0], 'd')
prank.getQuantiles(theproba, thequant)  # (3)
for i_q in range(nbquant.value):
    print("(theproba, thequant)[" + str(i_q) + "] = (" + str(theproba[i_q]) +
          ", " + str(thequant[i_q]) + ")")

allquant = ROOT.vector('double')()
prank.getQuantileVector(aproba, allquant)  # (4)
stdout.write("aproba = " + str(aproba.value) + ", allquant = ")
for quant_i in allquant:
    stdout.write(str(quant_i) + " ")
print("")

nbquant = 3
aproba = 0.8, aquant = 6.4
(theproba, thequant)[0] = (0.2, 1.6)
(theproba, thequant)[1] = (0.6, 4.8)
(theproba, thequant)[2] = (0.8, 6.4)
aproba = 0.8, allquant = 6.4 7.4 

The α-quantile can be evaluated by several ways:

Control variate

To estimate the α-quantile by control variate, you must use the computeQuantileCV method. The procedure to do this estimation is the following:

Importance sampling

To estimate the α-quantile by importance sampling, the method computeThreshold needs to be used. The procedure to make this estimation follows.

First, an object TImportanceSampling needs to be created. This object will allow the creation of a copy of the TDataServer where one of its attributes (sent in parameter) is replaced by a new attribute defined by its law (sent in parameter too).

tis = Sampler.TImportanceSampling(tds2, "rw", 
                                 DataServer.TNormalDistribution("rw_IS", 0.10, 0.015), nS)

And then, this new TDataServer must be collected via the getTDS method.

tdsis = tis.getTDS()	

A sampling needs to be generated for this new TDataServer:

sampis = Sampler.TSampling(tdsis,"lhs",nS)
sampis.generateSample()
tlfis = Launcher.TLauncherFunction(tdsis, "flowrateModel", "*", "Y_IS")
tlfis.setDrawProgressBar(False)
tlfis.run()	

Now, the probability of an output variable exceeding threshold can be computed with the computeThreshold method.

ISproba = tis.computeThreshold("Y_IS",seuil)	

For information, it is possible to compute the mean and standard deviation of this output variable.

ISmean = tis.computeMean("Y_IS")
ISstd = tis.computeStd("Y_IS")	

The Wilks quantile computation is an empirical estimation, based on order statistic which allows to get an estimation on the requested quantile, with a given confidence level , independently of the nature of the law, and most of the time, requesting less estimations than a classical estimation. Going back to the empirical way discussed in Section II.4.4.1: it consists, for a 95% quantile, in running 100 computations, ordering the obtained values and taking the one at either the 95-Th or 96-Th position (see the discussion on how to choose k in Section II.4.4.1). This can be repeated several times and will result in a distribution of all the obtained quantile values peaking at the theoretical value, with a standard deviation depending on the number of computations made. As it peaks on the theoretical value, 50% of the estimation are larger than the theoretical value while the other 50% are smaller.

In the following a quantile estimation of 95% will be considered with a requested confidence level of 95% (for more details on this method, see [metho]). If the sample does not exist yet, a possible solution is to estimate the minimum requested number of computations (which leads in our particular case to a sample of 59 events). Otherwise, one can ask Uranie the index of the quantile value for a given sample size, as such:

tds = DataServer.TDataServer("useless","foo")
quant=0.95
CL = 0.95
SampleSize=200
theindex=0
theindex = tds.computeIndexQuantileWilks(quant, CL, SampleSize)	  

The previous lines are purely informative, and not compulsory: the method implemented in Uranie to deal with the Wilks quantile estimation will start by calling these lines and complains if the minimum numbers of points is not available. In any case, the bigger the sample is, the more accurate the estimated value is. This value is finally determined using the method:

tds.addAttribute(DataServer.TNormalDistribution("attribute",0,1))
sampis = Sampler.TSampling(tds,"lhs", SampleSize)
sampis.generateSample()

from ctypes import c_double
value=c_double(0)
tds.estimateQuantile("attribute", quant, value, CL)	  

As stated previously, this is illustrated in a use-case macro which results in Figure XIV.6. There one can see results from two classical estimations of the 95% quantile. The distribution of their results is centered around the theoretical value. The bigger the sample is, the closer the average is to the theoretical value and the smaller the standard deviation is. But in any case, there is always 50% of estimation below and 50% above the theoretical value. Looking at the Wilks estimation, one can see that only 5% and 1% of the estimations are below the theoretical value respectively for the 95% and 99% confidence level distributions (at the price of smaller sample). With a larger sample, the standard deviation of the estimated value distribution for a 95% confidence level is getting smaller.

As for all methods above, this one has been modified so that it can handle constant-size vectors (at least given the pre-selection of event from the combination of the overall selection and the one provided in the method, as a second argument). As usual, the idea is to consider all elements of a vector independent from the other. If one considers the correlation matrix computed between two attributes, one being a scalar while the other one is a constant-sized vector with 10 elements, the resulting correlation matrix will be a 11 by 11 matrix.

Here are two examples of computeCorrelationMatrix calls, both using the tdstest.dat file already shown in Section II.3.4, which contains four attributes, three of which can be used here ( being a non constant-size vector, using it in this method will bring an exception error). In the following example, two correlation matrices are computed: the first one providing the correlation of both and attributes while the second focus on the former and only the second element of the latter.

"""
Example of correlation matrix computation for vector
"""
from rootlogon import DataServer

tdsop = DataServer.TDataServer("foo", "poet")
tdsop.fileDataRead("tdstest.dat")

# Consider a and x attributes (every element of the vector)
globalOne = tdsop.computeCorrelationMatrix("x:a")
globalOne.Print()

# Consider a and x attributes (cherry-picking a single element of the vector)
focusedOne = tdsop.computeCorrelationMatrix("x[1]:a")
focusedOne.Print()

This should lead to the following console return, where the first correlation matrix contains all pearson correlation coefficient (considering as a constant-size vector whose element are independent one to another) while the second on focus only on the second element of this vector (a vector's number start at 0). The following macro is shown in Section XIV.3.13.

4x4 matrix is as follows

     |      0    |      1    |      2    |      3    |
---------------------------------------------------------
   0 |          1      0.9449      0.6286       0.189 
   1 |     0.9449           1      0.8486         0.5 
   2 |     0.6286      0.8486           1      0.8825 
   3 |      0.189         0.5      0.8825           1 


2x2 matrix is as follows

     |      0    |      1    |
-------------------------------
   0 |          1         0.5 
   1 |        0.5           1