我用 R 编写了一个算法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}
来源1：https://github.com/wch/r-source/blob/trunk/src/nmath/qnorm.c#L101
来源2：https://github.com/wch/r-source/blob/843b6d3eb5c53ef585d3b0fe1bcb905e5c684de8/src/nmath/dpq.h
第一个参数（分位数p）的较小值可以R很好地近似输出。但是当我将其设置为 .7 时，它返回值nan。我不明白为什么。
70% 分位数可以像在 中一样正常计算R，但是这个表达式告诉我这个值不存在。
有人能发现我的错误吗？我已经考虑了好几天了。:-(
编辑我忘了提到我排除了使用log(p)值的选项，因此一些表达式，尤其是R_DT_qIv(p)和R_DT_CIv(p)被简化了。
答案1
对于任何感兴趣的人，代码最终都会产生预期的结果。代码如下。
\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}