You don't need PL/pgsql!


You don't need PL/pgsql to create functions that do some pretty sophisticated things. Here's a function written entirely in SQL that returns the inverse cumulative distribution function known in Microsoft Excel™ circles as NORMSINV. 
CREATE OR REPLACE FUNCTION normsinv(prob float8)
RETURNS float8
STRICT
LANGUAGE SQL
AS $$
WITH constants(a,b,c,d,p_low, p_high) AS (
    VALUES(
        ARRAY[-3.969683028665376e+01::float8 , 2.209460984245205e+02 , -2.759285104469687e+02 , 1.383577518672690e+02 , -3.066479806614716e+01 , 2.506628277459239e+00],
        ARRAY[-5.447609879822406e+01::float8 , 1.615858368580409e+02 , -1.556989798598866e+02 , 6.680131188771972e+01 , -1.328068155288572e+01],
        ARRAY[-7.784894002430293e-03::float8 , -3.223964580411365e-01 , -2.400758277161838e+00 , -2.549732539343734e+00 , 4.374664141464968e+00 , 2.938163982698783e+00],
        ARRAY[7.784695709041462e-03::float8 , 3.224671290700398e-01 , 2.445134137142996e+00 , 3.754408661907416e+00],
        0.02425::float8,
        (1 - 0.02425)::float8
    )
),
intermediate(p, q, r) AS (
    SELECT
        prob AS p,
        CASE
            WHEN prob < p_low AND prob > p_low THEN sqrt(-2*ln(prob))
            WHEN prob >= p_low AND prob <= p_high THEN prob - 0.5
            WHEN prob > p_high AND prob < 1 THEN sqrt(-2*ln(1-prob))
            ELSE NULL
        END AS q,
        CASE
            WHEN prob >= p_low OR prob <= p_high THEN (prob - 0.5)*(prob - 0.5)
            ELSE NULL
        END AS r
    FROM constants
)
SELECT
    CASE
        WHEN p  <  0 OR
             p  >  1 THEN 'NaN'::float8
        WHEN p  =  0 THEN '-Infinity'::float8
        WHEN p  =  1 THEN  'Infinity'::float8
        WHEN p  <  p_low THEN
            (((((c[1]*q+c[2])*q+c[3])*q+c[4])*q+c[5])*q+c[6]) /
             ((((d[1]*q+d[2])*q+d[3])*q+d[4])*q+1)
        WHEN p  >= p_low AND p <= p_high THEN
            (((((a[1]*r+a[2])*r+a[3])*r+a[4])*r+a[5])*r+a[6])*q /
            (((((b[1]*r+b[2])*r+b[3])*r+b[4])*r+b[5])*r+1)
        WHEN p  >  p_high THEN
            -(((((c[1]*q+c[2])*q+c[3])*q+c[4])*q+c[5])*q+c[6]) /
             ((((d[1]*q+d[2])*q+d[3])*q+d[4])*q+1)
        ELSE /* This should never happen */
            (p*0)/0 /* This should cause an error */
        END
FROM
    intermediate
CROSS JOIN
    constants
$$;

COMMENT ON FUNCTION normsinv(prob float8) IS $$This implementation is taken from https://stackedboxes.org/2017/05/01/acklams-normal-quantile-function/$$;
[Edit: There were some typos and wrong functions in the previous version]

Comments

  1. UnknownFebruary 23, 2019 at 10:28 AM

    Hi, it seems that your function returns very different result from Excel's NORMSINV function, for the same input. Can you please check it?

    ReplyDelete
  2. PLFebruary 24, 2019 at 9:50 PM

    Good example of the power of SQL. I think it can be marked as immutable.

    ReplyDelete
  3. Zlatko MatićMarch 8, 2019 at 4:24 AM

    Hi Deivid, "WHEN prob < p_low AND prob > p_low THEN sqrt(-2*ln(prob))" has to be changed to "WHEN prob >0 AND prob < p_low THEN sqrt(-2*ln(prob))"

    ReplyDelete
  4. JuniorAugust 8, 2019 at 4:11 AM


    In PG when I try to execute the function returns the error:
    Field "prob" does not exist.

    ReplyDelete

Post a Comment

Popular posts from this blog

Cool Stuff in PostgreSQL 10: Transition Table Triggers

Last time , we created a heirarchy of tables for logging, but we don't have anything they'd log, and we don't have a way to ensure the logging happens automagically.  This time, we'll fix those problems. First, a thing to log: CREATE TABLE foo(     id BIGSERIAL PRIMARY KEY,     t TEXT ); Next, let's make an efficient trigger function.  You may have seen one of these before, but it was probably row-based, so it couldn't be super efficient.  This one is statement-based, and uses the new transition table feature in PostgreSQL 10.  We only need to build one trigger function because our logging infrastructure is schema-agnostic.  We could have separated it into three trigger functions, but that's just more code to maintain.  It starts off as usual: CREATE OR REPLACE FUNCTION log() RETURNS TRIGGER LANGUAGE plpgsql AS $$ BEGIN     IF TG_OP NOT IN ('INSERT', 'UPDATE', 'DELETE') THEN    ...
Read more

Cool Stuff in PostgreSQL 10: Auto-logging

We started off by creating a logging infrastructure, then arranging for a single table to use it. Rather than repeat that work for each table, let's use a relatively although not completely new feature: `EVENT TRIGGER`. The idea here is that we fire a trigger on `CREATE TABLE` and see to it that the table is logged. We'll write the trigger first, even though in reality, we'd need to load the function it calls first. ```sql CREATE EVENT TRIGGER add_logger ON ddl_command_end WHEN tag IN ('create table') EXECUTE PROCEDURE add_logger(); COMMENT ON EVENT TRIGGER add_logger IS 'Ensure that each table which is not a log gets logged'; ``` The magic happens inside `add_logger()`, but it's magic we've already seen. First, we'll get the table's name and schema using `pg_event_trigger_ddl_commands()`, filtering out tables which are already log tables. The test here is crude and string-based, but we could easily go to schema-based...
Read more

Cool Stuff in PostgreSQL 10: Partitioned Audit Table

If you've dealt with PostgreSQL's answer to partitioning in the past, you may have put a lot of work into it.  Happily, you have a lot less to bother with as of PostgreSQL 10, which is not out yet as of this writing. Let's dive right in with a table that's a little like an audit trail. CREATE TABLE IF NOT EXISTS the_log (     "timestamp" timestamp with time zone DEFAULT now() NOT NULL,     "user" text NOT NULL DEFAULT CURRENT_USER,     action text NOT NULL,     table_schema text NOT NULL,     table_name text NOT NULL,     old_row jsonb,     new_row jsonb,     CONSTRAINT the_log_check CHECK (         (old_row IS NOT NULL)::integer +         (new_row IS NOT NULL)::integer > 0     ) ) PARTITION BY LIST(table_schema); COMMENT ON CONSTRAINT the_log_check ON the_log IS 'Make...
Read more