37. 数値計算・其の壱


R に用意されている「数値計算を行う道具」を紹介する.

数値演算誤差

R の内部では 10 進数を 2 進数に直して計算しているため,小数計算の種類によっては近似値(循環小数)になり,計算結果にゴミ(数値誤差)が生じる場合がある.このゴミを見えなくする場合は関数 round() を用いればよい.

 
 0.4-0.2-0.2                      # 普通の結果

[1] 0

 0.4-0.3-0.1                      # 結果にゴミが残る

[1] 2.775558e-17

 round( 0.4-0.3-0.1 )             # 関数 round() で丸める

[1] 0

丸めについて

R には実数を丸める以下の関数が用意されている.

関数

機能

round(x, digits = 0)

IEEE 式で丸めを行なう(デフォルトは digits = 0 で小数以下を丸める).

ceiling(x)

x 未満でない最小の整数を返す.

floor(x)

x 以上でない最大の整数を返す.

signif(x, digits = 6)

digits で指定された有効桁数(デフォルトは 6 桁)に丸める.

trunc(x)

x を「 0 へ向かって」整数化する.

関数 round() の丸め方は四捨五入とは微妙に異なる.IEEE では任意の桁位置における丸め処理法を次のように定めている( RjpWiki の記事,『工学のためのデータサイエンス入門』 間瀬 茂 他 著 (数理工学社) から引用した).

  1. 一番近い丸め結果候補が1つだけなら,その数に丸める.
  2. 一番近い丸め結果候補が2つある場合は,末尾が偶数のものに丸める(五入ばかりでなく五捨もあり得る!).
  3. 丸め処理は1段階で行なわなければならない.

この規則は丸めによる誤差が最小になる利点がある.ただし,規則 2. が適用される場合は通常の四捨五入とは異なる結果になる場合がある.また,規則 3. は,例えば 2.45 は直に 2 と丸めるべきであって,まず 2.5 としてから,次に 3 としてはならないことを主張している.

命令

round(2.4)

round(2.5)

round(3.5)

round(2.51)

round(3.51)

round(3.46)

結果

2

2

4

3

4

3

規則

1. : 四捨

2. : 五捨!

2. : 五入

1. : 五入

1. : 五入

1. と 3. : 四捨

ガンマ関数

以下のガンマ関数が用意されている.

関数

機能

beta(a, b)

ベータ関数 B(a,b) = (gamma(a) × gamma(b)) / gamma(a+b) .引数 a,b は零および負整数を除く数値(のベクトル).

lbeta(a, b)

ベータ関数の自然対数.直接計算しており,log(beta(a,b)) として計算しているわけではない.

gamma(x)

ガンマ関数.引数 x は零および負整数を除く数値(のベクトル).正の整数 n に対する階乗 n! は gamma(n+1) として計算する.

lgamma(x)

ガンマ関数の自然対数.直接計算しており,log(gamma(x)) として計算しているわけではない.

digamma(x)

lgamma の一階微分.

trigamma(x)

lgamma の二階微分.

tetragamma(x)

lgamma の三階微分.

pentagamma(x)

lgamma の四階微分.

choose(n, k)

二項係数(n 個から k 個を選ぶ場合の数).gamma(n+1) / ( gamma(k+1) × gamma(n-k+1) ) として計算しているので n, k が整数でなくても値は得られるが,結果は整数化される.

lchoose(n, k)

二項係数 choose(n,k) の自然対数.直接計算しており,log(choose(n,k)) として計算しているわけではない.

階乗の計算について,昔は n! を計算する場合は gamma(n+1) としたものだが,最近になって関数 factorial() というものが出てきたので,これを使って階乗計算をすることが出来るようになった.例えば 3! を計算する場合は以下のようにする.

 
 factorial(3)

[1] 6

prod() という関数を使っても計算できるが,0! = 1 とは計算してくれない.また,lgamma(n+1) と定義される関数 lfactorial(n) も用意された.

ベッセル関数

引数 x は非負値,次数 nu は実数(負の値の場合を含む)とする.

関数

機能

besselI(x, nu, expon.scaled=F)

変形された第 1 種のベッセル関数.

besselK(x, nu, expon.scaled=F)

変形された第 2 種のベッセル関数.

besselJ(x, nu)

第 1 種のベッセル関数.

besselY(x, nu)

第 2 種のベッセル関数.

gammaCody(x)

ガンマ関数.gamma() より計算は少し速いが gamma() より精度は落ちる.