Fonction d'estimateur de Theil-Sen en T-SQL

Réponses:

9

Je mentais quand j'ai dit que je ne pouvais pas le recoder en SQL. J'étais juste trop paresseux. Voici le code avec un exemple d'utilisation.

Le code est basé sur une bibliothèque Perl TheiSen , utilisant QuickMedian . Définissons un nouveau type de table pour transmettre facilement nos données à la procédure.

CREATE TYPE dbo.TheilSenInputDataTableType AS TABLE 
(
    ID INT IDENTITY(1,1),
    x REAL, 
    y REAL
)

Veuillez noter la colonne ID, qui est importante ici car notre solution utilise l'instruction CROSS APPLY afin d'obtenir une interprétation correcte de la boucle interne trouvée dans TheilSen.pm.

my ($x1,$x2,$y1,$y2);
foreach my $i(0 .. $n-2){
    $y1 = $y->[$i];
    $x1 = $x->[$i];
    foreach my $j($i+1 .. $n-1){
        $y2 = $y->[$j];
        $x2 = $x->[$j];

Nous aurons également besoin d'un nouveau type de données pour stocker un tableau de valeurs de type réel.

CREATE TYPE [dbo].[RealArray] AS TABLE(
    [val] [real] NULL
)

Voici la fonction f_QuickMedian , retournant la médiane pour un tableau donné. Le mérite en revient à Itzik Ben-Gan .

CREATE FUNCTION [dbo].[f_QuickMedian](@RealArray RealArray READONLY)
RETURNS REAL
AS
BEGIN
    DECLARE @Median REAL;
    DECLARE @QMedian REAL;

    SELECT @Median = AVG(1.0 * val)
    FROM
    (
        SELECT o.val, rn = ROW_NUMBER() OVER (ORDER BY o.val), c.c
        FROM @RealArray AS o
        CROSS JOIN (SELECT c = COUNT(*) FROM @RealArray) AS c
    ) AS x
    WHERE rn IN ((c + 1)/2, (c + 2)/2);

    SELECT TOP 1 @QMedian = val FROM @RealArray
    ORDER BY ABS(val - @Median) ASC, val DESC

    RETURN @QMedian
END

Et l' estimateur p_TheilSen :

CREATE PROCEDURE [dbo].[p_TheilSen](
      @TheilSenInput TheilSenInputDataTableType READONLY
    , @m Real OUTPUT
    , @c Real OUTPUT
)
AS
BEGIN
    DECLARE 
        @m_arr RealArray
      , @c_arr RealArray;       

    INSERT INTO @m_arr
        SELECT m
        FROM 
        (
            SELECT  
                t1.x as x1
                , t1.y as y1
                , t2o.x as x2
                , t2o.y as y2
                , t2o.y-t1.y as [y2 - y1]
                , t2o.x-t1.x as [x2 - x1]
                , CASE WHEN (t2o.x <> t1.x) THEN  CAST((t2o.y-t1.y) AS Real)/(t2o.x-t1.x) ELSE NULL END AS [($y2-$y1)/($x2-$x1)]
                , CASE WHEN t1.y = t2o.y THEN 0
                  ELSE
                    CASE WHEN t1.x = t2o.x THEN NULL
                        ELSE 
                        -- push @M, ($y2-$y1)/($x2-$x1);
                        CAST((t2o.y-t1.y) AS Real)/(t2o.x-t1.x)
                    END
                  END as m
            FROM @TheilSenInput t1
            CROSS APPLY
                    (
                    SELECT  t2.x, t2.y
                    FROM    @TheilSenInput t2
                    WHERE   t2.ID > t1.ID
                     ) t2o
        ) t
        WHERE m IS NOT NULL 

    SELECT @m = dbo.f_QuickMedian(@m_arr)

    INSERT INTO @c_arr
        SELECT y - (@m * x)
            FROM @TheilSenInput

    SELECT @c = dbo.f_QuickMedian(@c_arr)

END

Exemple:

DECLARE 
      @in TheilSenInputDataTableType
    , @m Real
    , @c Real

INSERT INTO @in(x,y) VALUES (10.79,118.99)
INSERT INTO @in(x,y) VALUES (10.8,120.76)
INSERT INTO @in(x,y) VALUES (10.86,122.71)
INSERT INTO @in(x,y) VALUES (10.93,125.48)
INSERT INTO @in(x,y) VALUES (10.99,127.31)
INSERT INTO @in(x,y) VALUES (10.96,130.06)
INSERT INTO @in(x,y) VALUES (10.98,132.41)
INSERT INTO @in(x,y) VALUES (11.03,135.89)
INSERT INTO @in(x,y) VALUES (11.08,139.02)
INSERT INTO @in(x,y) VALUES (11.1,140.25)
INSERT INTO @in(x,y) VALUES (11.19,145.61)
INSERT INTO @in(x,y) VALUES (11.25,153.45)
INSERT INTO @in(x,y) VALUES (11.4,158.03)
INSERT INTO @in(x,y) VALUES (11.61,162.72)
INSERT INTO @in(x,y) VALUES (11.69,167.67)
INSERT INTO @in(x,y) VALUES (11.91,172.86)
INSERT INTO @in(x,y) VALUES (12.07,177.52)
INSERT INTO @in(x,y) VALUES (12.32,182.09)


EXEC p_TheilSen @in, @m = @m OUTPUT, @c = @c OUTPUT

SELECT @m
SELECT @c

Retour:

m = 52.7079
c = -448.4853

Juste pour comparaison, la version perl renvoie les valeurs suivantes pour le même ensemble de données:

m = 52.7078651685394
c = -448.484943820225

J'utilise l'estimateur TheilSen pour calculer la métrique DaysToFill pour les systèmes de fichiers. Prendre plaisir!

kafe
la source
1

J'ai également vérifié, pour T-SQL, Oracle et les serveurs en général (trop complexes pour être écrits en SQL pur).

Cependant, vous pourriez être intéressé par ce (un ensemble scientifique / statistique pour Python). L'algorithme y est également implémenté et en Python. Python est un langage que les humains ont au moins une certaine chance de pouvoir comprendre, contrairement à Perl.

Votre question m'a intrigué et j'ai fouillé. Il existe des bibliothèques C et C ++ qui contiennent cet algorithme - et il est également disponible dans quelques packages R. Et le message de @srutzky semble également intéressant.

+1 pour une question intéressante BTW - et bienvenue sur le forum :-)

Vérace
la source