[pstricks] inverse standard normal

Matthias Ruess matthias.ruess at db.com
Thu Sep 23 11:59:57 CEST 2010

```Hi,

I need to draw the density function of the Vasicek distribution: This
density function has the following structure:

a * Exp[ (b * norminv(x) + c)^2 + d * norminv(x)^2  ] where a,b,c and d
are constants.

Since there exist an Exp-operator in pstricks, the only thing I need is a
norminv operator (an approximation of the inverse of the standard normal
distribution).
Is there someone who can help me to build this in pstricks?

On
http://home.online.no/~pjacklam/notes/invnorm/index.html#The_algorithm
you can find the algorithm in several programming languages.

Here a perl implementation:

sub ltqnorm (\$) {
#
# Lower tail quantile for standard normal distribution function.
#
# This function returns an approximation of the inverse cumulative
# standard normal distribution function.  I.e., given P, it returns
# an approximation to the X satisfying P = Pr{Z <= X} where Z is a
# random variable from the standard normal distribution.
#
# The algorithm uses a minimax approximation by rational functions
# and the result has a relative error whose absolute value is less
# than 1.15e-9.
#
# Author:      Peter John Acklam
# Time-stamp:  2000-07-19 18:26:14
# E-mail:      pjacklam at online.no
# WWW URL:     http://home.online.no/~pjacklam

my \$p = shift;
die "input argument must be in (0,1)\n" unless 0 < \$p && \$p < 1;

# Coefficients in rational approximations.
my @a = (-3.969683028665376e+01,  2.209460984245205e+02,
-2.759285104469687e+02,  1.383577518672690e+02,
-3.066479806614716e+01,  2.506628277459239e+00);
my @b = (-5.447609879822406e+01,  1.615858368580409e+02,
-1.556989798598866e+02,  6.680131188771972e+01,
-1.328068155288572e+01 );
my @c = (-7.784894002430293e-03, -3.223964580411365e-01,
-2.400758277161838e+00, -2.549732539343734e+00,
4.374664141464968e+00,  2.938163982698783e+00);
my @d = ( 7.784695709041462e-03,  3.224671290700398e-01,
2.445134137142996e+00,  3.754408661907416e+00);

# Define break-points.
my \$plow  = 0.02425;
my \$phigh = 1 - \$plow;

# Rational approximation for lower region:
if ( \$p < \$plow ) {
my \$q  = sqrt(-2*log(\$p));
return (((((\$c[0]*\$q+\$c[1])*\$q+\$c[2])*\$q+\$c[3])*\$q+\$c[4])*\$q+\$c[5])
/
((((\$d[0]*\$q+\$d[1])*\$q+\$d[2])*\$q+\$d[3])*\$q+1);
}

# Rational approximation for upper region:
if ( \$phigh < \$p ) {
my \$q  = sqrt(-2*log(1-\$p));
return
-(((((\$c[0]*\$q+\$c[1])*\$q+\$c[2])*\$q+\$c[3])*\$q+\$c[4])*\$q+\$c[5]) /
((((\$d[0]*\$q+\$d[1])*\$q+\$d[2])*\$q+\$d[3])*\$q+1);
}

# Rational approximation for central region:
my \$q = \$p - 0.5;
my \$r = \$q*\$q;
return (((((\$a[0]*\$r+\$a[1])*\$r+\$a[2])*\$r+\$a[3])*\$r+\$a[4])*\$r+\$a[5])*\$q
/
(((((\$b[0]*\$r+\$b[1])*\$r+\$b[2])*\$r+\$b[3])*\$r+\$b[4])*\$r+1);
}

Kind regards,
Matt

--

Informationen (einschließlich Pflichtangaben) zu einzelnen, innerhalb der EU tätigen Gesellschaften und Zweigniederlassungen des Konzerns Deutsche Bank finden Sie unter http://www.deutsche-bank.de/de/content/pflichtangaben.htm. Diese E-Mail enthält vertrauliche und/ oder rechtlich geschützte Informationen. Wenn Sie nicht der richtige Adressat sind oder diese E-Mail irrtümlich erhalten haben, informieren Sie bitte sofort den Absender und vernichten Sie diese E-Mail. Das unerlaubte Kopieren sowie die unbefugte Weitergabe dieser E-Mail ist nicht gestattet.

Please refer to http://www.db.com/en/content/eu_disclosures.htm for information (including mandatory corporate particulars) on selected Deutsche Bank branches and group companies registered or incorporated in the European Union. This e-mail may contain confidential and/or privileged information. If you are not the intended recipient (or have received this e-mail in error) please notify the sender immediately and delete this e-mail. Any unauthorized copying, disclosure or distribution of the material in this e-mail is strictly forbidden.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://tug.org/pipermail/pstricks/attachments/20100923/18917477/attachment.html>
```