原です。

間が空いてしまってすいません。

>高橋です.
>
>> あらかじめ表を作らずに逆関数求めるってのは、、、できなさそうだなあ。
>
>http://www5.airnet.ne.jp/tomy/cpro/sslib11.htm
>
>で近似できます.cですが,ただの数値計算なので書き直すのは簡単だと思いま
>す.精度は十分かと.

ありがとうございます。なるほど、関数が固定なので、とにかく計算できれば
なんでもありなんですねえ。

Ruby で書き直してみました。qnorm が -Infty..x の積分で、pnorm がその逆関
数です。

def qnorm(u)
  a = [1.24818987e-4, -1.075204047e-3, 5.198775019e-3,
       -0.019198292004, 0.059054035642, -0.151968751364,
       0.319152932694, -0.5319230073, 0.797884560593]
  b = [-4.5255659e-5, 1.5252929e-4, -1.9538132e-5,
       -6.76904986e-4, 1.390604284e-3, -7.9462082e-4,
       -2.034254874e-3, 6.549791214e-3, -0.010557625006,
       0.011630447319,-9.279453341e-3, 5.353579108e-3,
       -2.141268741e-3, 5.35310549e-4, 0.999936657524]
  
  u == 0.0 and return 0.5 
  y = u / 2.0
  y < -6.0 and return 0.0
  y > 6.0 and return 1.0
  y < 0.0 and y = - y
  if y < 1.0
    w = y * y
    z = a[0]
    (1...9).each do |i|
      z = z * w + a[i]
    end
    z *= (y * 2.0)
  else
    y -= 2.0
    z = b[0]
    (1...15).each do |i|
      z = z * y + b[i]
    end
  end

  u < 0.0 and return (1.0 - z) / 2.0
  return (1.0 + z) / 2.0
end

def pnorm(qn)
  b = [1.570796288, 0.03706987906, -0.8364353589e-3,
       -0.2250947176e-3, 0.6841218299e-5, 0.5824238515e-5,
       -0.104527497e-5, 0.8360937017e-7, -0.3231081277e-8,
       0.3657763036e-10, 0.6936233982e-12]

  if(qn < 0.0 || 1.0 < qn)
    $stderr.printf("Error : qn <= 0 or qn >= 1  in pnorm()!\n")
    return 0.0;
  end
  qn == 0.5 and return 0.0
  
  w1 = qn
  qn > 0.5 and w1 = 1.0 - w1
  w3 = -Math.log(4.0 * w1 * (1.0 - w1))
  w1 = b[0]
  1.upto 10 do |i|
    w1 += b[i] * w3**i;
  end
  qn > 0.5 and return Math.sqrt(w1 * w3)
  -sqrt(w1 * w3)
end


調べてみると十分な精度でした。実用としてはこれを使うのがいいですね。

もろに計算しているのを探したら、例の「C言語によるアルゴリズム事典」(技評)
に次のようなのがありました。正規分布の 0..z の値です。

SQ2PI = Math.sqrt(2 * Math::PI)

def norm_dist(z)
  z = z.to_f
  z2 = z * z
  t = q = z * Math.exp(-0.5 * z2) / SQ2PI
  3.step(199, 2) do |i|
    prev = q
    t *= z2 / i
    q += t
    if q == prev
      return q
    end
  end
  z > 0 ? 0.5 : -0.5
end

これは積分を部分積分で展開しているだけなので理解可能です(^^;。

その逆関数については載っていなかったので、これを使ってニュートン法
で求めたらどうなるかと、、、

def newton_a(y, ini, epsilon = 0.00000000001, limit = 30)
  x = ini
  limit.times do |i|
    prev = x
    f, df = yield(prev)
    x = (y - f)/df + prev
    break if (x - prev).abs < epsilon
  end
  x
end

def p_norm(y)
  newton_a(y, 0.0) {|x| [norm_dist(x), Math.exp(-0.5 * x**2) / SQ2PI]}
end

これでもまずまずのスピードで計算できました。


、、、ライブラリにまとめておいた方がいいなあ。しかしそれにしても
どこからどこまでの積分なのかすぐ分からなくなるのをなんとかしたい
ですねえ。

  normxXX_(z)  ... -Infty から z まで
  norm___x(z)  ... z から Infty まで
  norm__X_(z)  ... 0 から z まで

という記法はどうでしょう?

一応、正規分布表を載せます。__X_ の表です。

def mk_line(x, n)
    s = sprintf("%1.1f|", x)
    (0.00).step(0.09, 0.01) do |y|
      s << sprintf(" %1.#{n}f", yield(x + y))
    end
    s
end

def mk_tbl(s, e, n = 4, &b)
  printf("    " + (" "*(n-1) + "%1.2f")*10 + "\n", * (0..9).map {|x| x*0.01})
  (s.to_f).step(e.to_f, 0.1) do |x|
    puts mk_line(x, n, &b)
  end
end

alias norm__X_ norm_dist

mk_tbl(0.0, 3.1, 4) do |x|
  norm__X_(x)
end

       0.00   0.01   0.02   0.03   0.04   0.05   0.06   0.07   0.08   0.09
0.0| 0.0000 0.0040 0.0080 0.0120 0.0160 0.0199 0.0239 0.0279 0.0319 0.0359
0.1| 0.0398 0.0438 0.0478 0.0517 0.0557 0.0596 0.0636 0.0675 0.0714 0.0753
0.2| 0.0793 0.0832 0.0871 0.0910 0.0948 0.0987 0.1026 0.1064 0.1103 0.1141
0.3| 0.1179 0.1217 0.1255 0.1293 0.1331 0.1368 0.1406 0.1443 0.1480 0.1517
0.4| 0.1554 0.1591 0.1628 0.1664 0.1700 0.1736 0.1772 0.1808 0.1844 0.1879
0.5| 0.1915 0.1950 0.1985 0.2019 0.2054 0.2088 0.2123 0.2157 0.2190 0.2224
0.6| 0.2257 0.2291 0.2324 0.2357 0.2389 0.2422 0.2454 0.2486 0.2517 0.2549
0.7| 0.2580 0.2611 0.2642 0.2673 0.2704 0.2734 0.2764 0.2794 0.2823 0.2852
0.8| 0.2881 0.2910 0.2939 0.2967 0.2995 0.3023 0.3051 0.3078 0.3106 0.3133
0.9| 0.3159 0.3186 0.3212 0.3238 0.3264 0.3289 0.3315 0.3340 0.3365 0.3389
1.0| 0.3413 0.3438 0.3461 0.3485 0.3508 0.3531 0.3554 0.3577 0.3599 0.3621
1.1| 0.3643 0.3665 0.3686 0.3708 0.3729 0.3749 0.3770 0.3790 0.3810 0.3830
1.2| 0.3849 0.3869 0.3888 0.3907 0.3925 0.3944 0.3962 0.3980 0.3997 0.4015
1.3| 0.4032 0.4049 0.4066 0.4082 0.4099 0.4115 0.4131 0.4147 0.4162 0.4177
1.4| 0.4192 0.4207 0.4222 0.4236 0.4251 0.4265 0.4279 0.4292 0.4306 0.4319
1.5| 0.4332 0.4345 0.4357 0.4370 0.4382 0.4394 0.4406 0.4418 0.4429 0.4441
1.6| 0.4452 0.4463 0.4474 0.4484 0.4495 0.4505 0.4515 0.4525 0.4535 0.4545
1.7| 0.4554 0.4564 0.4573 0.4582 0.4591 0.4599 0.4608 0.4616 0.4625 0.4633
1.8| 0.4641 0.4649 0.4656 0.4664 0.4671 0.4678 0.4686 0.4693 0.4699 0.4706
1.9| 0.4713 0.4719 0.4726 0.4732 0.4738 0.4744 0.4750 0.4756 0.4761 0.4767
2.0| 0.4772 0.4778 0.4783 0.4788 0.4793 0.4798 0.4803 0.4808 0.4812 0.4817
2.1| 0.4821 0.4826 0.4830 0.4834 0.4838 0.4842 0.4846 0.4850 0.4854 0.4857
2.2| 0.4861 0.4864 0.4868 0.4871 0.4875 0.4878 0.4881 0.4884 0.4887 0.4890
2.3| 0.4893 0.4896 0.4898 0.4901 0.4904 0.4906 0.4909 0.4911 0.4913 0.4916
2.4| 0.4918 0.4920 0.4922 0.4925 0.4927 0.4929 0.4931 0.4932 0.4934 0.4936
2.5| 0.4938 0.4940 0.4941 0.4943 0.4945 0.4946 0.4948 0.4949 0.4951 0.4952
2.6| 0.4953 0.4955 0.4956 0.4957 0.4959 0.4960 0.4961 0.4962 0.4963 0.4964
2.7| 0.4965 0.4966 0.4967 0.4968 0.4969 0.4970 0.4971 0.4972 0.4973 0.4974
2.8| 0.4974 0.4975 0.4976 0.4977 0.4977 0.4978 0.4979 0.4979 0.4980 0.4981
2.9| 0.4981 0.4982 0.4982 0.4983 0.4984 0.4984 0.4985 0.4985 0.4986 0.4986
3.0| 0.4987 0.4987 0.4987 0.4988 0.4988 0.4989 0.4989 0.4989 0.4990 0.4990
3.1| 0.4990 0.4991 0.4991 0.4991 0.4992 0.4992 0.4992 0.4992 0.4993 0.4993