1

I coded an algorithm from R with pgfmath:

\documentclass{article}
\usepackage{tikz}
\usetikzlibrary{fpu}
\makeatletter
\pgfmathdeclarefunction{qnorm}{4}{%
  \begingroup
    \pgfkeys{/pgf/fpu,/pgf/fpu/output format=fixed}
    % #1 = p
    % #2 = sigma
    % #3 = mu
    % #4 = lower_tail (0 = true, 1 = false) % I'm not sure why it does this
    \pgfmathparse{#2 < 0 ? 1 : 0}       % if(sigma  < 0)    ML_ERR_return_NAN;
    \edef\sigma@negative{\pgfmathresult}    %
    %
    \pgfmathparse{#2 == 0 ? 1 : 0}      % if(sigma == 0)    return mu;
    \edef\sigma@zerocheck{\pgfmathresult}   %
    %
    \pgfmathparse{#4? exp(#1) : 1-exp(#1)}  % #define R_DT_qIv(p) lower_tail ? exp(p) : - expm1(p);
    \edef\p@under{\pgfmathresult}       % p_ = R_DT_qIv(p);
    %
    \pgfmathparse{\p@under - .5}        % q = p_ - 0.5;
    \edef\q{\pgfmathresult}         % 
    %
    \pgfmathparse{#4? -(1-exp(#1)) : exp(#1) }% #define R_DT_CIv(p) lower_tail ? -expm1(p) : exp(p);
    \edef\@tmp{\pgfmathresult}          % r = R_DT_CIv(p);
    %
    \pgfmathparse{abs(\q) <= .425 ? 1 : 0}  % if (fabs(q) <= .425) /* 0.075 <= p <= 0.925 */
    \edef\q@interval{\pgfmathresult}        % 
    %
    \pgfmathparse{
        \q@interval ?               % if (fabs(q) <= .425)
                .180625 - \q * \q       % r = .180625 - q * q;
                :
                ( \q > 0 ?              % if (q > 0)
                    sqrt(-ln(\@tmp))    % r = R_DT_CIv(p) & r = sqrt(-log(r))
                :           % else
                sqrt(-ln(\p@under)) % r = p_; & r = sqrt(-log(r))
            )
        }
    \edef\r{\pgfmathresult}
   % 
    \pgfmathparse{\r <= 5 ? \r - 1.6 : \r - 5}  % if (r <= 5.) (r += -1.6) else (r += -5)
    \edef\@@tmp{\pgfmathresult}
    \pgfmathparse{\@@tmp}           % add to r
    \edef\r{\pgfmathresult}
    \pgfmathparse{\q@interval ? 
    (#1 * (((((((\r * 2509.0809287301226727 +       % q is in the interval
                       33430.575583588128105) * \r + 67265.770927008700853) * \r +
                     45921.953931549871457) * \r + 13731.693765509461125) * \r +
                   1971.5909503065514427) * \r + 133.14166789178437745) * \r +
                 3.387132872796366608)
            / (((((((\r * 5226.495278852854561 +
                     28729.085735721942674) * \r + 39307.89580009271061) * \r +
                   21213.794301586595867) * \r + 5394.1960214247511077) * \r +
                 687.1870074920579083) * \r + 42.313330701600911252) * \r + 1))
                                    : ( \q > 0 ?        % q is positive but not within the interval
                                        ((((((((\r * 7.7454501427834140764e-4 + 
                       .0227238449892691845833) * \r + .24178072517745061177) *
                     \r + 1.27045825245236838258) * \r +
                    3.64784832476320460504) * \r + 5.7694972214606914055) *
                  \r + 4.6303378461565452959) * \r +
                 1.42343711074968357734)
                / (((((((\r *
                         1.05075007164441684324e-9 + 5.475938084995344946e-4) *
                        \r + .0151986665636164571966) * \r +
                       .14810397642748007459) * \r + .68976733498510000455) *
                     \r + 1.6763848301838038494) * \r +
                    2.05319162663775882187) * \r + 1))
                                        :           % q is not within the interval, nor is it positive
                                        ((((((((\r * 2.01033439929228813265e-7 +
                       2.71155556874348757815e-5) * \r +
                      .0012426609473880784386) * \r + .026532189526576123093) *
                    \r + .29656057182850489123) * \r +
                   1.7848265399172913358) * \r + 5.4637849111641143699) *
                 \r + 6.6579046435011037772)
                / (((((((\r *
                         2.04426310338993978564e-15 + 1.4215117583164458887e-7)*
                        \r + 1.8463183175100546818e-5) * \r +
                       7.868691311456132591e-4) * \r + .0148753612908506148525)
                     * \r + .13692988092273580531) * \r +
                    .59983220655588793769) * \r + 1.))
                    )
                 }
    \edef\val{\pgfmathresult}
    \pgfmathparse{\val}
    \edef\@@@tmp{\pgfmathresult}
    \pgfmathparse{(\sigma@negative ? sqrt(-1) : % if(sigma  < 0)    ML_ERR_return_NAN;
         (\sigma@zerocheck ? #3 :           % if(sigma == 0)    return mu;
            ( !(\q > 0) ?                   % if(q < 0.0) 
                #3+#2*-\@@@tmp :        % val = -val
                #3+#2*\@@@tmp           % return mu + sigma * val
            )
        )
        }
    \pgfmath@smuggleone\pgfmathresult
  \endgroup
}
\makeatother

\begin{document}
\pgfmathparse{qnorm(.05,2,2,1)}\pgfmathresult   % setting the first argument to .7 yields nan - I can't figure out why
\end{document}

Source 1: https://github.com/wch/r-source/blob/trunk/src/nmath/qnorm.c#L101

Source 2: https://github.com/wch/r-source/blob/843b6d3eb5c53ef585d3b0fe1bcb905e5c684de8/src/nmath/dpq.h

Small values for the first argument (quantile p) approximate the R output nicely. But when I set this to .7, it returns value nan. I can't figure out why.

The 70% quantile can be calculated as normally in R, but this expression tells me this value doesn't exist.

Can anyone spot my mistake ? I've sat on this for days. :-(

EDIT I forgot to mention that I exluded the option to use log(p) values, so some expressions, notably R_DT_qIv(p) and R_DT_CIv(p) are simplified.

1 Answer 1

1

For anyone interested, the code finally produces the intended result. Code's below.

\documentclass{article}
\usepackage{tikz}
\usetikzlibrary{fpu}
\makeatletter
\pgfmathdeclarefunction{qnorm}{4}{%
  \begingroup
    \pgfkeys{/pgf/fpu,/pgf/fpu/output format=fixed}
    \pgfmathparse{#1}
    \edef\p{\pgfmathresult}
    \pgfmathparse{#2}
    \edef\sigma{\pgfmathresult}
    \pgfmathparse{#3}
    \edef\mu{\pgfmathresult}
    %
    \pgfmathparse{\sigma < 0 ? 1 : 0}           % if(sigma  < 0)    ML_ERR_return_NAN;
    \edef\sigma@neg{\pgfmathresult}
    \pgfmathparse{\sigma == 0 ? 1 : 0}          % if(sigma == 0)    return mu;
    \edef\sigma@zero{\pgfmathresult}
    %
    \pgfmathparse{\p - .5}                  % q = p_ - 0.5;
    \edef\q{\pgfmathresult}
    %
    \pgfmathparse{abs(\q) <= .425 ? 1 : 0}      % if (fabs(q) <= .425) /* 0.075 <= p <= 0.925 */
    \edef\q@i{\pgfmathresult}               % 
    %
    \pgfmathparse{                      % q is in interval [.075 ; .925]
      \q@i ?                        % if {(fabs(q) <= .425)
        .180625 - \q * \q               % r = .180625 - q * q;}
      :                     % /* r = min(p, 1-p) < 0.075 */
        sqrt(-ln(min(\p,(1-\p))))                       % 
    }
    \edef\r@{\pgfmathresult}
   % 
    \pgfmathparse{\q@i ? 
      \r@ % r = .180625 - q * q;
    : 
      (\r@ <= 5 ?  % if (r <= 5.)
        \r@ - 1.6 %
      : 
        \r@ - 5
      )}
    \edef\r{\pgfmathresult}
    \pgfmathparse{\q@i ? 
    (\q * (((((((\r * 2509.0809287301226727 +
                       33430.575583588128105) * \r + 67265.770927008700853) * \r +
                     45921.953931549871457) * \r + 13731.693765509461125) * \r +
                   1971.5909503065514427) * \r + 133.14166789178437745) * \r +
                 3.387132872796366608)
            / (((((((\r * 5226.495278852854561 +
                     28729.085735721942674) * \r + 39307.89580009271061) * \r +
                   21213.794301586595867) * \r + 5394.1960214247511077) * \r +
                 687.1870074920579083) * \r + 42.313330701600911252) * \r + 1.))
                                    : ( \r@ <= 5 ?  % if (r <= 5.)      %
                                        ((((((((\r * 7.7454501427834140764e-4 +     %
                       .0227238449892691845833) * \r + .24178072517745061177) *     %
                     \r + 1.27045825245236838258) * \r +                        %
                    3.64784832476320460504) * \r + 5.7694972214606914055) *     %
                  \r + 4.6303378461565452959) * \r +                        %
                 1.42343711074968357734)                                %
                / (((((((\r *                                       %
                         1.05075007164441684324e-9 + 5.475938084995344946e-4) *     %
                        \r + .0151986665636164571966) * \r +                    %
                       .14810397642748007459) * \r + .68976733498510000455) *       %
                     \r + 1.6763848301838038494) * \r +                     %
                    2.05319162663775882187) * \r + 1))                      %
                                        :       % /* very close to  0 or 1 */
                                        ((((((((\r * 2.01033439929228813265e-7 +
                       2.71155556874348757815e-5) * \r +
                      .0012426609473880784386) * \r + .026532189526576123093) *
                    \r + .29656057182850489123) * \r +
                   1.7848265399172913358) * \r + 5.4637849111641143699) *
                 \r + 6.6579046435011037772)
                / (((((((\r *
                         2.04426310338993978564e-15 + 1.4215117583164458887e-7)*
                        \r + 1.8463183175100546818e-5) * \r +
                       7.868691311456132591e-4) * \r + .0148753612908506148525)
                     * \r + .13692988092273580531) * \r +
                    .59983220655588793769) * \r + 1.))
                    )
                 }
    \edef\val{\pgfmathresult}
    \pgfmathparse{\q@i ? \val :
    (\q > 0 ? \val : -\val)}                    % if(q < 0.0)  --> val = -val ( in diff. r case)
    \edef\@tmp{\pgfmathresult}
    \pgfmathparse{\p > 1 ? sqrt(-1) :
      (\p < 0 ? sqrt(-1) :
        \sigma@neg ? sqrt(-1) :             % if(sigma  < 0)    ML_ERR_return_NAN;
         (\sigma@zero ? \mu :               % if(sigma == 0)    return mu;
           (#4 ? \mu+\sigma*\@tmp :         
                     \mu-\sigma*\@tmp)              % return mu + sigma * val
         )
       )
    }
    \pgfmath@smuggleone\pgfmathresult
  \endgroup
}
\makeatother
\begin{document}
\pgfmathparse{qnorm(.4,2,0,0)}\pgfmathresult    %
\end{document}
3
  • what did you do to fix it?
    – Ben Bolker
    Jun 10, 2014 at 20:45
  • @BenBolker There was an issue with the natural logarithm. \p@under evaluated to negative values and that caused issues. I changed around definitions of the min/max functions to plainly reflect the desired value rather than taking versatility for expansion of the function with log(p) expressions into account. After that I shuffled around a few expressions (mostly about the last \pgfmathparse{}) by determining each step in both R and in my algorithm to make them match. The lower/upper tail problem I kind of solved with a workaround that works way better in LaTeX than in C.
    – 1010011010
    Jun 10, 2014 at 20:51
  • Namely, by mirroring \val based on the given expression for #4.
    – 1010011010
    Jun 10, 2014 at 20:56

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .