有时候不一定所有的数学算法都有类库，碰到一些棘手的也只能硬上了。希望能给大家一些帮助，以下是我的实现
理论为一篇论文，详细说明如何进行johnson变换
大体过程为1.设集合Z={0.25，,026,,...,1.25} 针对每个z,取-3z,-z,z,3z 对应标准正态分布的累积分布函数值q1，q2，q3,q4
2.计算q1 q2 q3 q4 在源数据中的四分位数(percentile) 为x1 x2 x3 x4 ,计算QR = (x4 - x3)*(x2 - x1)/(Math.Pow(x3 - x2, 2));
3.计算
double lamda, ita, eplsion, gama;
double n = x2 - x1;
double p = x3 - x2;
double m = x4 - x3;
double q0 = .5*(x2 + x3);
double mp = m / p;
double np = n / p;
4.针对QR的大小分情况讨论
QR>1 源数据使用Su变换
ita = 2*(double)z/(Fn.ArCosh(.5*(mp+np)));
gama = ita*Fn.ArSinh((np-mp)/(2*Math.Sqrt(mp*np-1)));
lamda = 2*p*Math.Sqrt(mp*np-1)/((mp+np - 2)*Math.Sqrt(mp+np + 2));
eplsion = q0+p*(np-mp)/(2*(mp+np-2));
（公式 gama + ita*Fn.ArSinh((x - epsilon)/lamda))
然后使用AD检验 得出p值
QR<1 源数据使用SB变换
double pm = p/m;
double pn = p/n;
ita = (double)z / Fn.ArCosh(.5 * Math.Sqrt(((1 + pm) * (1 + pn))));
gama = ita* Fn.ArSinh(((pn-pm)*Math.Sqrt((1+pm)*(1+pn)-4) / (2*(pm*pn-1))));
lamda = p * Math.Sqrt((((1 + pm) * (1 + pn) - 2) * ((1 + pm) * (1 + pn) - 2) - 4)) / (pm * pn - 1);
eplsion = q0 - .5*lamda + p*(pn-pm) / (2*(pm*pn-1));
SB变换公式： gama + ita*Math.Log(x - epsilon);
然后使用AD检验 得出p值
QR=1 源数据使用Sl变换
ita = 2*(double)z/(Math.Log(mp));
gama = ita*Math.Log(Math.Abs(mp - 1)/(p*Math.Sqrt(mp)));
eplsion = q0 - .5*p*(mp + 1)/(mp - 1);
变换公式： gama + ita*Math.Log(x - epsilon);
然后使用AD检验 得出p值
6.自此我们计算所有Z的对应变换的p值，取p值最大的那个Z所对应的变换就是我们要的johnson变换
看起来算法还是蛮复杂的，里面有些函数可以使用类库，一般也比较常见。
小心，经过验算结果还是和minitab软件结果基本一致的。