суббота, 11 августа 2018 г.

Статистические справочные таблицы на калькуляторе

Задача. Есть калькулятор, но нет под рукой статистических таблиц. Например, нужны таблицы критических точек распределения Стьюдента для вычисления доверительного интервала. Взять компьютер с Excel? Не спортивно.
Большая точность не нужна, можно воспользоваться приближенными формулами. Идея приведённых ниже формул состоит в том, что преобразованием аргумента все распределения можно так или иначе свести к нормальному. Аппроксимации должны обеспечивать как вычисление кумулятивной функции распределения, так и расчет обратной к ней функции.
Начнём с нормального распределения. $$\Phi(z)=P=\frac{1}{2}\left[1+\mathrm{erf}\left(\frac{z}{\sqrt{2}}\right)\right]$$ $$z=\Phi^{-1}(P)=\sqrt{2}\cdot\mathrm{erf}^{-1}(2P-1)$$
Для него требуется вычислить функцию \(\mathrm{erf}(x)\) и обратную к ней. Я воспользовался приближением [1]: $$ \mathrm{erf}(x)=\mathrm{sign}(x)\cdot\sqrt{1-\exp\left(-x^{2}\cdot\frac{\frac{4}{\pi}+ax^{2}}{1+ax^{2}}\right)} $$ $$ \mathrm{erf}^{-1}(x)=\mathrm{sign}(x)\cdot\sqrt{-t_{2}+\sqrt{t_{2}^{2}-\frac{1}{a}\cdot\ln t_{1}}} $$ где \(t_1\) и \(t_2\) вспомогательные переменные: $$ t_{1}=1-x^{2},\:t_{2}=\frac{2}{\pi a}+\frac{\ln t_{1}}{2} $$ а константа \(a=0.147\). Ниже дан код на языке Octave.
function y = erfa(x)
  a  = 0.147;
  x2 = x**2; t = x2*(4/pi + a*x2)/(1 + a*x2);
  y  = sign(x)*sqrt(1 - exp(-t));
endfunction
function y = erfinva(x)
  a  = 0.147; 
  t1 = 1 - x**2; t2 = 2/pi/a + log(t1)/2;
  y  = sign(x)*sqrt(-t2 + sqrt(t2**2 - log(t1)/a));
endfunction

function y = normcdfa(x)
  y = 1/2*(1 + erfa(x/sqrt(2)));
endfunction
function y = norminva(x)
  y = sqrt(2)*erfinva(2*x - 1);
endfunction
Теперь, когда есть функции нормального распределения, приведём аргумент и вычислим t-распределение Стьюдента [2]: $$ F_{t}(x,n)=\Phi\left(\sqrt{\frac{1}{t_{1}}\cdot\ln(1+\frac{x^{2}}{n})}\right) $$ $$ t=F_{t}^{-1}(P,n)=\sqrt{n\cdot\exp\left(\Phi^{-1}(P)^{2}\cdot t_{1}\right)-n} $$ где вспомогательная переменная \(t_1\) есть $$ t_{1}=\frac{n-1.5}{(n-1)^{2}} $$
function y = tcdfa(x,n)
  t1 = (n - 1.5)/(n - 1)**2;
 y = normcdfa(sqrt(1/t1*log(1 + x**2/n)));
endfunction
function y = tinva(x,n)
  t1 = (n - 1.5)/(n - 1)**2;
  y  = sqrt(n*exp(t1*norminva(x)**2) - n);
endfunction
Идея приближенного вычисления распределения \(\chi^{2}\) наглядно представлена формулами [3]: $$ \sigma^{2}=\frac{2}{9n},\:\mu=1-\sigma^{2} $$ $$ F_{\chi^{2}}(x,n)=\Phi\left(\frac{\left(\frac{x}{n}\right)^{1/3}-\mu}{\sigma}\right) $$ $$ \chi^{2}=F_{\chi^{2}}^{-1}(P,n)=n\cdot\left(\Phi^{-1}(P)\cdot\sigma+\mu\right)^{3} $$
function y = chi2cdfa(x,n)
  s2 = 2/9/n; mu = 1 - s2;
  y  = normcdfa(((x/n)**(1/3) - mu)/sqrt(s2));
endfunction
function y = chi2inva(x,n)
 s2 = 2/9/n; mu = 1 - s2;
  y = n*(norminva(x)*sqrt(s2) + mu)**3;
endfunction
Распределение Фишера (для \(n/k\geq3\) и \(n\geq3\)) находится в два шага. Сначала аргумент преобразуется к вычислению распределения Фишера через распределение \(\chi^{2}\) [4], а его мы уже знаем, как вычислить. $$ \sigma^{2}=\frac{2}{9n},\:\mu=1-\sigma^{2} $$ $$ \lambda=\frac{2n+k\cdot x/3+(k-2)}{2n+4k\cdot x/3} $$ $$ F_{f}(x;k,n)=\Phi\left(\frac{\left(\lambda\cdot x\right)^{1/3}-\mu}{\sigma}\right) $$ Найдём обратную функцию, решив квадратное уравнение. $$ q=\left(\Phi^{-1}(P)\cdot\sigma+\mu\right)^{3} $$ $$ b=2n+k-2-4/3\cdot kq $$ $$ D=b^{2}+8/3\cdot knq $$ $$ x=F_{f}^{-1}(P;k,n)=\frac{-b+\sqrt{D}}{2k/3} $$
function y = fcdfa(x,k,n)
  mu = 1-2/9/k; s = sqrt(2/9/k);
  lambda = (2*n + k*x/3 + k-2)/(2*n + 4*k*x/3);
  normcdfa(((lambda*x)**(1/3)-mu)/s)
endfunction
function y = finva(x,k,n)
  mu = 1-2/9/k; s = sqrt(2/9/k);
  q = (norminva(x)*s + mu)**3;
  b = 2*n + k-2 -4/3*k*q;
  d = b**2 + 8/3*k*n*q;
  y = (sqrt(d) - b)/(2*k/3);
endfunction

Список литературы

  1. Sergei Winitzki. A handy approximation for the error function and its inverse. February 6, 2008.
  2.  Gleason J.R. A note on a proposed Student t approximation // Computational statistics & data analysis. – 2000. – Vol. 34. – №. 1. – Pp. 63-66.
  3.  Wilson E.B., Hilferty M.M. The distribution of chi-square // Proceedings of the National Academy of Sciences. – 1931. – Vol. 17. – №. 12. – Pp. 684-688.
  4.  Li B. and Martin E.B. An approximation to the F-distribution using the chi-square distribution. Computational statistics & data analysis. – 2002. Vol. 40. – №. 1. pp. 21-26.

Комментариев нет:

Отправить комментарий