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 6 (of 6 total)
You must be logged in to reply to this topic. Login to reply
This website stores cookies on your computer.
These cookies are used to improve your website experience and provide more personalized services to you, both on this website and through other media.
To find out more about the cookies we use, see our Privacy Policy