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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.