inverse normal

[Andre]

Hi everyone,

I am using excel formula like "NORMSINV(0.875)". Anyone know convert to TCL function?

Thanks in advance


Cheers!

[Arjen Markus]

The math::statistics package does contain a procedure that returns the inverse of the error function. It is used (or was used) internally only, so it is not documented or exported, but you can use it nonetheless. See Inverse-cdf-normal in this package (source file: pdf_stat.tcl)

Regards,

Arjen

[Andre]

I'm trying to verify the calculation for the Inverse Cumulative Standard
Normal Distribution Function presented by the Excel function NORMSINV(p).

For instance: The Excel function =NORMSINV(0.9626) returns 1.78.

Could you show me a manual calculation for this? I would like to do in TCL scripting.

Thanks in advance

Cheers!

[Christian Gollwitzer]

Am 25.09.14 08:38, schrieb Andre:

I've also had this problem once and translated the algorithm from

http://home.online.no/~pjacklam/notes/invnorm/#Pseudo_code_for_rational_approximation

into Tcl:

> #################################################
> ## ::tcl::mathfunc::invnorm
> #
> # inverse normal distribution in expr namespace
> # Uses rational approximation from
> # http://home.online.no/~pjacklam/notes/invnorm/#Pseudo_code_for_rational_approximation
> # relative precision 1.2*10^-9 in the full range
> #
> ###################################################
>
> proc ::tcl::mathfunc::invnorm {p} {
> # inverse normal distribution
> # rational approximation from
> # http://home.online.no/~pjacklam/notes/invnorm/#Pseudo_code_for_rational_approximation
> # precision 1.2*10^-9
>
> if {$p<=0 || $p>=1} {
> error "Domain error (invnorm)"
> }
> # Coefficients in rational approximations.
> 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
>
> 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
>
> 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
>
> set d1 7.784695709041462e-03
> set d2 3.224671290700398e-01
> set d3 2.445134137142996e+00
> set d4 3.754408661907416e+00
>
> # Define break-points.
>
> set p_low 0.02425
> set p_high [expr {1-$p_low}]
>
> # Rational approximation for lower region.
>
> if {$p < $p_low} {
> set q [expr {sqrt(-2*log($p))}]
> set x [expr {((((($c1*$q+$c2)*$q+$c3)*$q+$c4)*$q+$c5)*$q+$c6) / \
> (((($d1*$q+$d2)*$q+$d3)*$q+$d4)*$q+1)}]
> return $x
> }
>
> # Rational approximation for central region.
>
> if {$p <= $p_high} {
> set q [expr {$p - 0.5}]
> set r [expr {$q*$q}]
> set x [expr {((((($a1*$r+$a2)*$r+$a3)*$r+$a4)*$r+$a5)*$r+$a6)*$q / \
> ((((($b1*$r+$b2)*$r+$b3)*$r+$b4)*$r+$b5)*$r+1)}]
> return $x
> }
>
> # Rational approximation for upper region.
>
>
> set q [expr {sqrt(-2*log(1-$p))}]
> set x [expr {-((((($c1*$q+$c2)*$q+$c3)*$q+$c4)*$q+$c5)*$q+$c6) /
> (((($d1*$q+$d2)*$q+$d3)*$q+$d4)*$q+1)}]
> return $x
>
> }
>
> puts [expr {invnorm(0.9626)}]

Running this, I get back 1.7816887672191446

Christian

[Andre]

Hi Christian,

Very fantastic Chris but I think that still having incorrect this issue. I was try change the data set below;

set a1 2.515517
set a2 0.802853
set a3 0.010328

set b1 1.432788
set b2 0.189269
set b3 0.001308

Please advice.

Cheers!

[Christian Gollwitzer]

Hi Andre,

please do contextual quoting like this (
http://en.wikipedia.org/wiki/Posting_style#Interleaved_style ). It is
much easier to read.

Am 25.09.14 12:55, schrieb Andre:

> On Thursday, September 25, 2014 2:14:44 PM UTC+7, Christian Gollwitzer wrote:
>>> puts [expr {invnorm(0.9626)}]
>>
>> Running this, I get back 1.7816887672191446
>

> Very fantastic Chris but I think that still having incorrect this issue. I was try change the data set below;
>
> set a1 2.515517
> set a2 0.802853
> set a3 0.010328
>
> set b1 1.432788
> set b2 0.189269
> set b3 0.001308

I don't understand. What have you tried? Have you, by any chance,
changed the numbers inside the invnorm function? Those are values
internal to the algorithm, you should not change them at all.

After you source the script I posted just as is, you can compute any
inverse normal distribution by doing

expr {invnorm($p)}

where p is set to the probability.

Christian

[Andre]

Hi Christian,

Everything is going well. Thank you very much, Crist!

[Andre]

Hi Christian,

Do you know the expression for "critical values of Z from grubb's test? You could check it out http://graphpad.com/quickcalcs/grubbs1/.

Thanks in advance

Cheers!

[Christian Gollwitzer]

Hi Andre,

Am 29.09.14 07:57, schrieb Andre:

> Do you know the expression for "critical values of Z from grubb's
> test? You could check it out
> http://graphpad.com/quickcalcs/grubbs1/.

sorry, never heard about Grubb's test. Since this is a new question,
don't post it as a response under your old question. Open a new thread
instead. And since this is not a Tcl related question, but rather a
statistics one, please ask in sci.math.

Christian

[Andre]

Copied, Christian!

I will go there to find it. Many thanks for this, Christ!