作者別: ARQ

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

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

• 日本人全体のうち, あるもの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番目の手続きに戻る.

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

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

2018年11月16日にキログラムの定義を変更することが CGPM (国際度量衡総会)で可決され, 本日2019年5月20日により全世界的に施行されることとなった.

いままでのキログラムの定義は次のように物質に依存する定義であった。

• 水1リットルの重量
• 大気圧下で氷が溶けつつある温度における水1リットルの重量
• 最大密度(液温摂氏4度)における蒸留水1立方デシメートル(1リットル)の質量
• アルシーブ原器 (有体に言うと白金塊) の質量
• IPK (国際キログラム原器: 白金とイリジウムの合金)の質量

水はいろんなものを自然に溶かし込むし、金属も徐々にいろんなものを吸着して変質するので不変で安定したものというにははばかられた.

それでもIPKは誤差1億分の5という精度ではあったが, 近年の科学技術の発達により, これでも精度は足りないという状況であった.

今回の改訂は, 特殊相対性理論から得られた式 E=mc² と光子エネルギーの式 E=hν から m=hν/c² を得るので, 今までは実験値であったプランク定数を真の意味で定数とすることにより キログラムが定義されるということになる。

ただ, この定義に至る裏側で, 日本やドイツを中心とする科学者たちは (IPKと比べれば)手軽に作成できる原器を作り, その精度を裏付ける役割を担った. その原器はケイ素(シリコン)結晶である. このあたりの話は産業技術総合研究所のウェブサイトが詳しい.

非常に大きな数を扱うための概念として対数を高校数学で学習する. 例えば, 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 の範囲にあるといえる.

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

今年もよろしくお願いいたします.

昨年の積み残しがあるので今年は何とかクリアしたいです. 何とか次のものは達成したいところ.

• 数学検定1級合格
• 統計検定1級合格
• TOEICで人並みの点数を取れるようになるか英検準1級をとれるように頑張る
• 5kmを40分以内に走り切れる体力を取り戻す

学習の成果はこのブログにもちょっとずつ載せていきますので, 小難しい記事が増えるかもしれません...

今後ともよろしくお願いいたします.

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

ユークリッドの互除法とは, 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) とすることができる.

具体例