カテゴリー: 数学

ベイズの定理とは事後確率と事前確率の変換を表したものである. ここでは最も簡単なモデルを取り扱ってみよう.

ここでは以下の仮定をする.

• 日本人全体のうち, あるものMをもっているという比率を π とする.
• あるものをもっている日本人に対し質問Qを行ったときに期待した回答(私は確かにMをもっている)が得られる確率をπ_Aとする.
• あるものをもっていない日本人に対し質問Qを行ったときに期待した回答(私は確かにMをもっていない)が得られる確率をπ_Sとする.

これらの仮定のもとである日本人に対し質問Qを行ったところ, 私は確かにMをもっているという回答を得た. その人が本当に M をもっている確率を求めてみよう.

結論から言うと, これはベイズの定理によって計算できる.

実際に M をもっているという事象を H, 質問Qにより私は確かにMをもっているという回答を得るという事象をAとする.

求める確率は Pr(H|A) であり, これはベイズの定理により Pr(A|H)Pr(H)/(Pr(A|H)Pr(H)+Pr(A|￢H)Pr(￢H)) に等しい.

これに仮定で設定した値を代入すると (π_Aπ)/(π_Aπ+(1-π_S)(1-π)) を得る.

最近話題になっている例の奴について当てはめると...

某所にて, (要約すると)ArcTan(5/12) および ArcTan(3/4) を度数法で表せ(小数以下切り捨て) なる問題を出題した. そのときは円周率 π が 333/106<π<355/113 と近似できるというヒントを付けた.

これと, sin(x)<x<tan(x) (第1象限の角に限る) という大小関係をうまく用いてもらい, 適当な評価をしてもらう算段でいた.

実際に, sin(x)=S, tan(x)=T は正確に求められる状況であり, π₁<π<π₂ という条件を与えれば, x [rad]= θ [deg] という関係においては 180S/π₂<θ<180T/π₁ という関係が成り立つ.

ただ, この評価方法は0°<θ<15° 程度でないと大雑把すぎるので, 4倍角や5倍角の正弦や正接を計算させる必要があった。

ところが, 必要であれば2倍角の計算をすればいいようなもっといい近似方法があったようだ. 屈折率の式でおなじみのヴィレブロルト・スネルが見出した式である, 3sin(x)/(2+cos(x))<x<(2sin(x)+tan(x))/3 である (以降, この不等式の左辺をL, 右辺をUとする).

これにより、角度の評価は 180L/π₂<θ<180U/π₁ となる.

何といっても, その誤差は x⁵に比例するということが大きいし係数も小さい.

下限の比較では sin(x)=x-x³/6+Ο(x⁵) に対し, 3sin(x)/(2+cos(x))=x-x⁵/180+Ο(x⁷) となり, 上限の比較では tan(x)=x+x³/3+Ο(x⁵) に対し, (2sin(x)+tan(x))/3=x+x⁵/20+Ο(x⁷) となる。

この方法で当初の問題を解いてみると...

θ=ArcTan(5/12) のとき, 15/38<ArcTan(5/12)<185/468 は度数法ではこの時点で15900/703<θ<20905/923 と評価できる. この時点で小数第1位まで確定する. (θ=22.6...°)

θ=ArcTan(3/4) のとき, 9/14<θ<13/20 は度数法ではこの時点で9540/259<θ<13221/355 と評価できる. ざっと, θ=37°±0.1° といったところか.

さらに2倍角の公式で得た角と直角との差を評価すると 21/74<π/2-2Arctan(3/4)<511/1800 となり, 122767/3330<θ<50475/1369 を得る. なんと, これで小数第2位まで確定する (θ=36.86...°)

9回目の受験にして漸く実用数学技能検定1級に合格することができました!

初めて挑戦したのは6年前. 直前に受験した準1級を満点で合格したのだから1級も何とかなるだろうと思ってチャレンジして, その結果が惨敗だったときの悔しさを胸に, 続けて2回受験するもまるで歯が立たず...

転居で環境が変わったところで気軽に受験できなくなった(個人受験会場が近くになかった)ために諦めていたところ, 次の転居ではまた個人受験会場が近くなったことから再起を決意.

昨年から受験をしたもののやはり合格には届かなかったが, 解法が全く思い浮かばないということがなくなったので、もう少しで手が届くという感触は得ていた。

それでも, 前回の結果が逆ボーダー （1次が4.5/7 (合格ライン 5/7), 2次が (2.3/4 (合格ライン 2.5/4)) だったときは心が折れかけた...

が, ここまでできているのだから何とかなると気合を入れ直したところ, 今度こそ合格ラインを超えることができたのである.

記事の都合上, 解釈が面倒な数式を書くので, 同じ式を表現した WolframAlpha へのリンクを設定しています.

連分数とは、分数の分母にさらに分数が含まれるような数値表現である.

例えば, 1+2/(3+4/(5+6/7))というようなものである.

特に、分子の部分がすべて 1 になっているような連分数のことを正則連分数といい, 単に連分数といった時には正則連分数であることを暗黙で示していることがある.

正則連分数について, a₀+1/(a₁+1/(a₂+1/a₃)) を [a₀; a₁, a₂, a₃] と表現する. 各 a_i は正の整数 (a₀ は0や負の値を含めた整数) である.

ある実数を(正則)連分数に変換することを連分数展開というが, 有理数ならば有限項の連分数で, 無理数ならば無限項の連分数で表現されることがわかっている.例えば, 355/113=[3; 7, 16], π=[3; 7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, 14, 2, 1, 1, 2, 2, 2, 2, 1, 84, 2, 1, 1, 15, 3, 13, 1, ...] となる.

無理数を分数(整数比)で近似する場合, 連分数展開を途中で打ち切って真分数や仮分数に戻すという手段が使われる。 そうして現れる分子や分母の値は, それ以下の整数を用いたあらゆる整数比(分数)よりもよりよい近似となることが知られている.

実際に, 例に挙げた 355/113=[3; 7, 16]=[3; 7, 15, 1] は 355 以下の整数を用いた分数の中で最良の円周率の近似である.

さて, 正則連分数展開は次の手順で行われる.

1. 正則連分数展開を行いたい数を ω₀ とし, n を 0 とする.
2. ω_n の整数部分を a_n とする.
3. ω_n の小数部分が 0 に等しい場合は, a_k (k≥n+1) は存在しない, すなわち, 連分数表現は a_n の項で打ち切る.
4. ω_n+1 を ω_n=a_n+1/ω_n+1 を満たすように定める.
5. n に 1 を加えて2番目の手続きに戻る.

この作業を行う数式処理ソフトををいくつか紹介する.

連分数近似を実装している数式処理ソフトの紹介

非常に大きな数を扱うための概念として対数を高校数学で学習する. 例えば, 2¹⁰ は 10³ とほぼ等しい. 10 の累乗に値を換算したときの指数に相当する数を, その元の数の常用対数 をとった値というのである. これを式で表したものが log₁₀2¹⁰≈3 という関係である.

なお, 3dB(decibel( で 2 倍という関係はこの対数計算から得られるのである.

これを一般化すると, 1 ではない正の実数 a, 任意の実数 r について対数(logarithm(を次のように定義するということである。

a^r=R ⇔ log_aR=r

a=10 のときは常用対数, a=e (e: ネイピア数) のときは自然対数, a=2 のときは二進対数と特別に呼ばれる.

log_aR という表記であるが, 文脈から明らかな場合は底 a はしばしば省略される.

ここでよく話題になる話として, 様々な整数の常用対数をとるといくらになるのか? というものがある. いくつかの整数に対して, その概算値を出してみよう.

log(1) は正確に 0 である. なぜならば, すべての 0 ではない実数 R について R⁰=1 が成り立つからである.

log(2) は冒頭で話題にした 2¹⁰≈10³ の関係から 10log(2)≈3 ⇔ log(2)≈0.3 となる. 厳密には 2¹⁰>10³ であるから log(2)>0.3 である. 上からの評価もしたい. 1000<1024<1250 の関係を使うと. 常用対数をとると 3<10log(2)<=4-3log(2) となるから log(2)<4/13<0.3077 が得られるが、 ちょっとこれではおおざっぱすぎる. ここは天下り的ではあるが, 2⁹³<10²⁸ という関係を使う. すると, log(2)<28/93<0.3011 を得る.

log(3) についてもここでは近似値を求めることになる. 80<81 より 1+3log(2)<4log(3) から log(3)>19/40=0.475 という評価ができる. また, 20000>19683=3⁹ ⇔ 4+log(2)>9log(3) から log(3)<4/9+log(2)/9=400/837<0.4779

これらの値を用いて log(4)=2log(2), log(5)=1-log(2), log(6)=log(2)+log(3), log(8)=3log(2), log(9)=2log(3) が概算できる.

1桁の整数の残り log(7) は様々な近似が検討されている. 例えば, 48<49<50 ⇔ 4log(2)+log(3)<2log(7)<2-log(2) より 67/80=0.8375<log(7)<17/20=0.85 という簡易的な評価がある. もうちょっと精度を上げるならば 7⁴=2401>2400 ⇔ 4log(7)>2+3log(2)+log(3) より log(7)>27/32=0.84375, 48×1024/1000>49 ⇔ 14log(2)+log(3)-3>2log(7) より log(7)<1417/1674<0.8465 さらには, 10⁶<2²×3⁶6×7³ と 2¹⁰×7⁸<3¹⁰×10⁵ を評価すれば小数第3位まで確定できる.

最後に log_e(10) を評価する. そうすると, 自然対数と常用対数の変換に必要な定数が得られるのである. ただ, こちらについては積分計算により求めるのが一般的だろう. ∫₁²(1/x)dx=log_e(2) という関係があるからだ. 被積分関数は積分区間内において凸関数なので、台形近似では収束がちょっと遅い. シンプソンの公式なら劇的に早く収束する. 具体的には (1/(6n))(1+2Σ_j=1^n-11/(1+2j/(2n))+4Σ_j=1ⁿ(1/(1+(2j-1)/(2n)))+1/2) これを必要な精度まで計算して log₁₀(2) で割ると log_e(10) が得られる. シンプソンの公式は4次収束なので、n=2 でも3桁の精度は出る. log_e(2)≈0.693 なので log_e(10) は 2.25 から 2.31 の範囲にあるといえる.

おまけ: 連分数による近似

世界最古に発見されたアルゴリズムのひとつとして知られるユークリッドの互除法は, 歴史上の価値のみならず, 公約数という素因数が絡む計算を素因数分解せずに効率的に求めることができるという点で非常に有用な手法となっています.

ユークリッドの互除法とは, 2つの0ではない整数a, b に対して, 次の手続を繰り返すことで a と b の最大公約数を求めることができます.

1. b=0 ならば a が求めるべき最大公約数である. アルゴリズムはここで終了する.
2. c≡a (mod b) を計算し, a に b を代入, b に c を代入する.
3. 最初に戻る.

ラメの定理により, アルゴリズムの繰り返し回数は a, b のうち小さい方の 10 進桁数を d とすれば, 高々 5d 回となることが知られており, 5d 回程度の繰り返しが発生する例はフィボナッチ数列(最初の2項を1としたもの)の隣接2項を対象とした場合が知られています.

もしも, これを素因数分解で解決しようとすると, a, b のうち小さい方の正の平方根の数だけ「試しに割り算して割り切れるかを確かめる」作業が必要となります。

つまり, a≥bのもとで割算の回数を比較すると, ユークリッドの互除法ならば 2.2log(b), 素因数分解するなら √b と見積もることができます. b が大きい場合, どちらの方が効率がいいのかは明らかです.

この方法を拡張すると不定方程式 ax+by=gcd(a,b) の整数解を導出することができます. これが拡張されたユークリッドの互除法です. 拡張されたユークリッドの互除法の手順は次のとおりです.

1. r₀=a, s₀=1, t₀=0, r₁=b, s₁=0, t₀=1, i=1 とする.
2. r_i+1=r_i-1-q_ir_i を 0≤r_i+1≤abs(r_i), q_i は整数となるように定める.
3. s_i+1=s_i-1-q_is_i, t_i+1=t_i-1-q_it_i とする.
4. r_i+1=0 の場合は x=s_i, y=t_i であるので計算を終了する. r_i+1≠0 の場合は次の項目に進む.
5. i に 1 を加えて 2. の項目へ戻る.

ユークリッドの互除法の使用例

個人的メモ

法 m における整数 a のモジュラ逆数 a^-1 (mod m) とは, aa^-1≡1 (mod m) を満たす a^-1≡x (mod m) なる整数 x の合同類である.

法 m における整数 a のモジュラ逆数が存在するための必要十分条件は gcd(m, a)=1 である.

a^-1 (mod m) を求める最も汎用的な方法は拡張されたユークリッドの互除法を用いることである. つまり, a, m が所与であるもとでの ax-qm=gcd(a,m)=1 を満たす x が a^-1 (mod m) である.

他にも, オイラーのトーシェント関数を用いて a^-1≡a^φ(m)-1 (mod m) としたり, カーマイケル関数を用いて a^-1≡a^λ(m)-1 (mod m) と求める方法があるが, これは m を因数分解する必要があるので効率的な導出方法ではない. これらの特殊形として, m が素数 p であった場合は, a^-1≡a^p-2 (mod p) とすることができる.

具体例

個人的なメモ

• フェルマーの小定理とは次の関係を言う: 素数 p と p とは互いに素な整数 a について a^p-1≡1 (mod p) が成り立つ.
• 数論におけるオイラーの定理とは次の関係を言う: 任意の正の整数 n と n とは互いに素な整数 a について a^φ(n)≡1 (mod n) が成り立つ.
  • φ(n) はオイラーのトーシェント関数であり、これは n 以下の n とは互いに素な正の整数の数となる。
  • n=Π_k=1^mp_k^b_k と素因数分解できたならば, φ(n)=nΠ_k=1^m(1-1/p_k)
• カーマイケルの定理とは次の関係を言う: φ(n) は a^m≡1(mod n) を満たす最小の m であるとは限らない.
  最小の m はカーマイケル関数λ(n)で与えられる.
  • カーマイケル関数は次のように再帰的に定義される.
    • λ(2^e)=1 (e=1), 2 (e=2), 2^e-2 (e≧3)
    • 奇素数 p について λ(p^e)=p^e-1(p-1)
    • n=Π_k=1^mp_k^b_k と素因数分解できたならば λ(n)=lcm(λ(p₁^b₁),...,λ(p_m^b_m)). ただし, lcm は引数たちの最小公倍数(least common multiple)を表す。

  n≡1 (mod λ(n)) を満たすような n をカーマイケル数という.

以上によれば, 2ⁿ≡1 (mod 15) を満たす n は, フェルマーの小定理では導出できない (n=3×5)が,
オイラーの定理によれば n=8 のときに成立することが分かり, カーマイケルの定理によれば, 最小の整数 n は 4 である.

三角形(頂点を A, B, C とし、その対辺の長さをそれぞれ a, b, c とする)の面積を求める方法には、次のものが知られている。

• (底辺)×(高さ)÷2
• (1/2)ab sin∠C
• ヘロンの公式(Heron's formula): √(s(s-a)(s-b)(s-c)), ただし, s=(a+b+c)/2

では、これを四角形に拡張できないものかと考えていたが、すでにブレートシュナイダーの公式(Bretschneider's formula)が存在することが知られている(ブラーマグプタの公式(Brahmagupta's formula)では四角形が円に内接するという条件が必要であるためプログラムには適さない)。四角形(頂点を順にA, B, C, D とし, AB=p, BC=q, CD=r, DA=s)に対応するブレートシュナイダーの公式は次のとおりである。

√((T-p)(T-q)(T-r)(T-s)-pqrscos²((A+C)/2))), ただし, T=(p+q+r+s)/2

ただし、今回の場合は頂点の座標から p, q, r, s, ∠A, ∠C を導出する必要がある.よって, 与えられた頂点の座標を四角形の頂点の順番に並べ替えなければならない.

いろんな方法を模索している中で、次の記述を発見した。

1. x座標が小さい順番に並べ、それぞれの点を P₁, P₂, P₃, P₄ とする.
• P_1y<P_2yのとき
  • P_3y>P_3yとなるように
    P₃, P₄を入れ替える.
• P_1y>P_2yのとき
  • P_3y<P_3yとなるように
    P₃, P₄を入れ替える.

一見してよさそうに見えるが, ∠P₃が平角よりも大きくなる場合はうまくいかない. それどころか,
P₂とP₃を入れ替えてもなお四角形となる場合があることに気づいてしまった.
なので, 座標だけを与えられても凹四角形の場合はうまくいかないことが判明した...

よって、頂点の座標から四角形の面積を計算する場合は、次のような設計を考えなければならない.

• 入力者は頂点の座標を四角形の頂点の順序となるように入力する.
• 面積を算出する側は与えられた頂点座標から四角形が描けることを確認してから角度と辺の長さを計算しブレートシュナイダーの公式を適用する.

ということで、現在、面積を算出する側は与えられた頂点座標から四角形が描けることを確認する手段を思索中である.
おそらく, 内角の和が 2π になることは四角形であることの必要十分条件であることをいえばいいはずなのだが...

もうちょっと考えてみた

XREAサーバーが増強されたことを受けて加筆しています. 不要になった部分はあえて見え消しにしています. テキストをそのままコピーすると余計な部分まで取り込む虞があります.

車輪の再発明... というよりも, プログラム作成練習ということで, 統計でよく利用される分布の統計数値表や統計量などを計算するプログラムを作成してみた. 知る人ぞ知る青木繁伸先生が公開されているプログラムのクローンである. 使用したプログラム言語は C++ で, Boost ライブラリを利用した. Boost の開発者が発表した特殊関数計算にかかる近似誤差は, 青木先生が公開している AWK プログラムを直接 C++ に移植したものよりも小さいので, より正確な値が出力されるものと期待できる.

現在間借りしている XREA は, CGI プログラムとしてマシンネイティブなバイナリを利用することが認められている. ただし, 設置者の環境で行ったものをアップロードする必要がある. サーバーの環境を調べてみると, Linux 32ビットLinux 64ビット環境で動いていることが分かったので...

1. Windwos 用のプログラムバイナリを使うことができない. Linux 用のプログラムを生成する必要がある
2. 筆者の環境は64ビットアーキテクチャであるが, 32ビット版のバイナリを生成する必要がある.64ビット版のバイナリか32ビット版のバイナリを用意する必要がある.

という問題が発生した.

どうやって対応したか