February 1, 2005 at 11:21 am
Since SQL does not seem to have the 'inverse of the cummulative standarized normal distribution' function (NORMSINV), I am trying to find a way to do this.
Does anyone have a clue?
Thanks in advance!
-Steve
February 2, 2005 at 6:00 am
Not my work but I found this on the web:
http://home.online.no/~pjacklam/notes/invnorm/
Good luck,
Rick
create or replace procedure stdnormal_inv(p in float, u out float)
is type doubleArray is varray(6) of float;
a doubleArray := doubleArray( -3.969683028665376e+01, 2.209460984245205e+02, -2.759285104469687e+02, 1.383577518672690e+02, -3.066479806614716e+01, 2.506628277459239e+00 );
b doubleArray := doubleArray( -5.447609879822406e+01, 1.615858368580409e+02, -1.556989798598866e+02, 6.680131188771972e+01, -1.328068155288572e+01 );
c doubleArray := doubleArray( -7.784894002430293e-03, -3.223964580411365e-01, -2.400758277161838e+00, -2.549732539343734e+00, 4.374664141464968e+00, 2.938163982698783e+00 );
d doubleArray :=doubleArray( 7.784695709041462e-03, 3.224671290700398e-01, 2.445134137142996e+00, 3.754408661907416e+00 );
q float; t float;
begin if (p = 0.0) then u:=-10000;
return;
end if;
if (p = 1.0) then u:=10000;
return;
end if;
q := least(p,1-p);
if (q > 0.02425) then /* Rational approximation for central region. */ u := q-0.5;
t := u*u;
u := u*(((((a(1)*t+a(2))*t+a(3))*t+a(4))*t+a(5))*t+a(6)) /(((((b(1)*t+b(2))*t+b(3))*t+b(4))*t+b(5))*t+1);
else /* Rational approximation for both tail regions. */ t := sqrt(-2*ln(q));
u := (((((c(1)*t+c(2))*t+c(3))*t+c(4))*t+c(5))*t+c(6)) /((((d(1)*t+d(2))*t+d(3))*t+d(4))*t+1);
end if;
/* The relative error of the approximation has absolute value less than 1.15e-9. One iteration of Halley's rational method (third order) gives full machine precision... */ stdnormal_cdf(u,t);
t:=t-q;
/* error */ t := t*sqrt(2*3.141592654)*exp(u*u/2);
/* f(u)/df(u) */ u := u-t/(1+u*t/2);
/* Halley's method */ if(p>0.5) then u:=-u;
end if;
return;
end stdnormal_inv;/
February 2, 2005 at 11:55 am
Rick, Thanks. I need to run this on a MS SQL server, so I did some translating, but can't get past a few things.
Maybe you or someone else can help me out. I will post my modified (MS SQL) version of the code. I left most of the original comments in so that I could understand it better.
Look at the section I have commeted with /* This is the part I don't understand. What is stdnormal_cdf and what variable gets set? Is there a SQL equivalent? */
create function stdnormal_inv(@p float)
RETURNS float
as
BEGIN
declare @a1 float
declare @a2 float
declare @a3 float
declare @a4 float
declare @a5 float
declare @a6 float
set @a1 = -3.969683028665376e+01
set @a2 = 2.209460984245205e+02
set @a3 = -2.759285104469687e+02
set @a4 = 1.383577518672690e+02
set @a5 = -3.066479806614716e+01
set @a6 = 2.506628277459239e+00
declare @b1 float
declare @b2 float
declare @b3 float
declare @b4 float
declare @b5 float
set @b1 = -5.447609879822406e+01
set @b2 = 1.615858368580409e+02
set @b3 = -1.556989798598866e+02
set @b4 = 6.680131188771972e+01
set @b5 = -1.328068155288572e+01
declare @c1 float
declare @c2 float
declare @c3 float
declare @c4 float
declare @c5 float
declare @c6 float
set @c1 = -7.784894002430293e-03
set @c2 = -3.223964580411365e-01
set @c3 = -2.400758277161838e+00
set @c4 = -2.549732539343734e+00
set @c5 = 4.374664141464968e+00
set @c6 = 2.938163982698783e+00
declare @d1 float
declare @d2 float
declare @d3 float
declare @d4 float
set @d1 = 7.784695709041462e-03
set @d2 = 3.224671290700398e-01
set @d3 = 2.445134137142996e+00
set @d4 = 3.754408661907416e+00
declare @q float
declare @t float
declare @U float
if (@p = 0.0)
begin
set @U=-10000
return @U
end
else if (@p = 1.0)
begin
set @U=10000
return @U
end
else
begin
/*set @q = least(@p,1-@p) Can't do "least" in SQL,
had to do something different */
if @p > (1-@p)
begin
set @q = (1-@q)
end
else
begin
set @q = @p
end
if (@q > 0.02425)
begin
/* Rational approximation for central region. */
set @U = @q-0.5
set @t = @U*@u
set @U = @U*(((((@a1*@t+@a2)*@t+@a3)*@t+@a4)*@t+@a5)*@t+@a6)/(((((@b1*@t+@b2)*@t+@b3)*@t+@b4)*@t+@b5)*@t+1)
end
else
begin
/* Rational approximation for both tail regions. */
set @t = sqrt(-2*log(@q))
set @U = (((((@c1*@t+@c2)*@t+@c3)*@t+@c4)*@t+@c5)*@t+@c6)/((((@d1*@t+@d2)*@t+@d3)*@t+@d4)*@t+1)
/* The relative error of the approximation has absolute value less
than 1.15e-9. One iteration of Halley's rational method (third
order) gives full machine precision... */
/* This is the part I don't understand. Is there a SQL equivalent?
What is stdnormal_cdf and what variable gets set? */
stdnormal_cdf(@u,@t)
set @t=@t-@q /* error */
set @t = @t*sqrt(2*3.141592654)*exp(@u*@u/2) /* f(u)/df(u) */
set @U = @u-@t/(1+@u*@t/2) /* Halley's method */
end
if(@p>0.5)
begin
set @U=-@u
end
end
return @U
END
June 28, 2006 at 4:24 am
Steve
It looks to me like if you follow the link in the previous post to http://home.online.no/~pjacklam/notes/invnorm/ you'll see two functions -
the one you have translated (well done !) and another one ...
stdnormal_cdf
This also needs translating, but I believe it would set @t in your code.
Good Luck
Dom
September 26, 2006 at 11:35 am
Thank you for NORMSINV function ! It's work !!
Here is a SQL code for SYBASE 12.5 version:
IF EXISTS
(SELECT * FROM sysobjects WHERE type = 'P' AND name = 'GET_NORMSINV')
BEGIN
PRINT 'Dropping Procedure GET_NORMSINV'
DROP Procedure GET_NORMSINV
END
GO
'Creating Procedure GET_NORMSINV'
GO
CREATE Procedure
GET_NORMSINV
@p
float,
@result
float OUTPUT
AS
/******************************************************************************
** File:
** Name: GET_NORMSINV
** Desc: '-- Equivalent of Excel''s NORMSINV function
** The relative error of the approximation has absolute value less
** than 1.15e-9. One iteration of Halley's rational method (third
** order) gives full machine precision...
**
** This template can be customized:
**
** Return values:
**
** Called by:
**
** Parameters:
**
**
*******************************************************************************
** Change History
*******************************************************************************
** Date: Author: Description:
** -------- -------- -------------------------------------------
**
*******************************************************************************/
declare @a1 float
declare @a2 float
declare @a3 float
declare @a4 float
declare @a5 float
declare @a6 float
declare @b1 float
declare @b2 float
declare @b3 float
declare @b4 float
declare @b5 float
declare @c1 float
declare @c2 float
declare @c3 float
declare @c4 float
declare @c5 float
declare @c6 float
declare @d1 float
declare @d2 float
declare @d3 float
declare @d4 float
declare @plow float
declare @phigh float
declare @q float
declare @r float
select @a1 = -39.6968302866538
select @a2 = 220.946098424521
select @a3 = -275.928510446969
select @a4 = 138.357751867269
select @a5 = -30.6647980661472
select @a6 = 2.50662827745924
select @b1 = -54.4760987982241
select @b2 = 161.585836858041
select @b3 = -155.698979859887
select @b4 = 66.8013118877197
select @b5 = -13.2806815528857
select @c1 = -7.78489400243029E-03
select @c2 = -0.322396458041136
select @c3 = -2.40075827716184
select @c4 = -2.54973253934373
select @c5 = 4.37466414146497
select @c6 = 2.93816398269878
select @d1 = 7.78469570904146E-03
select @d2 = 0.32246712907004
select @d3 = 2.445134137143
select @d4 = 3.75440866190742
select @plow=.02425
select @phigh=1-@plow
if (@p<@plow)
begin
select @q = Sqrt(-2 * Log(@p))
select @result=(((((@c1 * @q + @c2) * @q + @c3) * @q + @c4) * @q + @c5) * @q + @c6) / ((((@d1 * @q + @d2) * @q + @d3) * @q + @d4) * @q + 1)
end
else
begin
if (@p<@phigh)
begin
select @q =@p - 0.5
select @r = @q * @q
select @result= (((((@a1 * @r + @a2) * @r + @a3) * @r + @a4) * @r + @a5) * @r + @a6) * @q / (((((@b1 * @r + @b2) * @r + @b3) * @r + @b4) * @r + @b5) * @r + 1)
end
else
begin
select @q = Sqrt(-2 * Log(1 - @p))
select @result= -(((((@c1 * @q + @c2) * @q + @c3) * @q + @c4) * @q + @c5) * @q + @c6) / ((((@d1 * @q + @d2) * @q + @d3) * @q + @d4) * @q + 1)
end
end
GO
GRANT EXEC ON
GET_NORMSINV TO PUBLIC
GO
September 27, 2006 at 3:01 am
Here is an equivalent of Excel's NORMSDIST function (wrote for SYbase 12.5):
IF EXISTS (SELECT * FROM sysobjects WHERE type = 'P' AND name = 'GET_NORMSDIST')
BEGIN
PRINT 'Dropping Procedure GET_NORMSDIST'
DROP Procedure GET_NORMSDIST
END
GO
PRINT 'Creating Procedure GET_NORMSDIST'
GO
CREATE Procedure GET_NORMSDIST
@x float,
@result float OUTPUT
AS
/******************************************************************************
** File:
** Name: GET_NORMSDIST
** Desc: '-- Equivalent of Excel''s NORMSDIST function
** The relative error of the approximation has absolute value less
** than 1.15e-120
** This template can be customized:
**
** Return values:
**
** Called by:
**
** Parameters:
**
**
*******************************************************************************
** Change History
*******************************************************************************
** Date: Author: Description:
** -------- -------- -------------------------------------------
**
*******************************************************************************/
declare @L float
declare @k float
declare @dCND float
declare @pi float
declare @a1 float
declare @a2 float
declare @a3 float
declare @a4 float
declare @a5 float
select @L = 0.0
select @k = 0.0
select @dCND = 0.0
select @a1 = 0.31938153
select @a2 = -0.356563782
select @a3 = 1.781477937
select @a4 = -1.821255978
select @a5 = 1.330274429
select @pi = 3.1415926535897932384626433832795
select @L = Abs(@x)
if @L >= 30
begin
if sign(@x) = 1
select @result = 1
else
select @result = 0
end
else
begin
-- perform calculation
select @k = 1.0 / (1.0 + 0.2316419 * @L)
select @dCND = 1.0 - 1.0 / Sqrt(2 * @pi) * Exp(-@L * @L / 2.0) *
(@a1 * @k + @a2 * @k * @k + @a3 * POWER(@K, 3.0) + @a4 * POWER(@K, 4.0) + @a5 * POWER (@K, 5.0))
if (@x < 0)
select @result = 1.0 - @dCND
else
select @result = @dCND
end
GO
GRANT EXEC ON GET_NORMSDIST TO PUBLIC
GO
Viewing 6 posts - 1 through 5 (of 5 total)
You must be logged in to reply to this topic. Login to reply