- 追加された行はこの色です。
- 削除された行はこの色です。
#contents
*演習の目的 [#b76d5bdb]
演習１の目的は，forループを活用して，各種統計値を算出することである．堆積学では特に角度データに関する統計値を多用するが，それらはExcelの組み込み関数などでは計算できない．したがって，自分でプログラミングを行い，角度データ統計値を算出するツールを作り出す必要がある．
*参考書 [#z26d2210]
通常の統計学に関しては無数の日本語の参考書があるので，ここで挙げることはしない．一方，角度データの統計学に関しては日本語の良い教科書は存在しない（砕屑物の研究法には多少の記述がある）．英語では，
|Graham Borradaile, 2003, Statistics of Earth Science Data - Their Distribution in Time, Space, and Orientation. Springer.|
が非常に優れている．他には，
|John C. Davis, 2002, Statistics and Data Analysis in Geology Third Edition. John Wiley & Sons.|
|Graham J.G. Upton and Bernard Fingleton, 1989, Spatial Data Analysis by Example. Volume 2 Categorical and Directional Data. John Wiley & Sons.|
が役に立つ．
*スカラーデータの平均値と標準偏差 [#rdad27ac]
まず，基本的な記述統計値である「平均値」を計算しよう．スカラーデータの平均値\(\bar{d}\) は，
\[
\bar{d} = \frac{ \sum {d_i} } {n}
\]
で定義される．ここで，\(d_i\) は\(i\) 番目のデータ，\(n\) はデータの個数を意味する．
:問題1|
平均値を求めるプログラムを書いてみよう．計算のためのデータは&ref(data1.txt);を利用すること．
:問題1の解答例|
Scilabを使っている場合は&ref(calc_average.sce);を参照．Javaの場合は&ref(StatsMain.java);を参照．
次に，スカラーデータの標準偏差を計算する．標準偏差\(s\) は次の式で定義される：
\[
s = \sqrt{ \frac{\sum{ (d_i - \bar{d})^2 } }{n} }
\]
標準偏差が大きいほどデータセット数値のばらつきが大きい．一方，標本標準偏差（不偏標準偏差）と呼ばれる統計値もある．これは，データセットの母集団の標準偏差に対する不偏推定量となっている．標本標準偏差\(\sigma\) は以下の式で求められる：
\[
\sigma = \sqrt{ \frac{\sum{ (d_i - \bar{d})^2 } }{n - 1} }
\]
標準偏差と標本標準偏差は\(n\) で割るか\(n-1\) で割るかの違いであり，\(n\)が十分に大きい場合はほとんど値に違いがなくなる．単に標準偏差といった場合は標本標準偏差を指していることもあるので注意が必要である．
:問題2|
標準偏差を求めるプログラムを書いてみよう．計算のためのデータは前回同様&ref(data1.txt);を利用すること．
:問題2の解答例|
Scilabを使っている場合は&ref(calc_stdev.sce);を参照．
*スカラーデータ平均値の95%信頼区間 [#gb26cce6]
非常に大きい数の母集団から\(n\) 個の標本を取り出しているとする．このとき，標本から得られる平均値\(\bar{d}\) は母集団の真の平均値\(\mu\) とは必ずしも一致しない．この場合，以下の数式から得られる\(t\) の確率密度分布は自由度\(n-1\) の\(t\) 分布に従うことが知られている：
\[
\label{t_def}
t = \frac{\bar{d} - \mu}{ \sigma / \sqrt{n} }
\]
ここで\(t\) 分布の具体的な求め方を示すことはしない（興味がある人は統計学の教科書を参照のこと）．ただし，式\(\ref{t_def}\)を変形すると，以下の式が得られることは重要である：
\[
\bar{d} - \mu = t \cdot { \frac{\sigma} {\sqrt{n}} }.
\]
すなわち，標本から得られる平均値と‟真の”平均値のズレの大きさの確率密度分布は，\(t\) の確率密度分布に\(\sigma / \sqrt{n}\)を掛けたものに等しいのである．（この\(s_e=\sigma / \sqrt{n}\) を標準誤差と呼ぶ）．つまり，\(1 - \alpha\)%の確率で真の（母集団）の平均値が含まれている範囲は，
\[
\label{t_confint}
\bar{d} \pm t_{\alpha/2,n-1} \cdot s_e
\]
ということになる．ここで，\(t_{\alpha/2,n-1}\) は，\(\alpha/2\)% の面積を含む\(t\) 分布の両側境界値である．
さて，一見求めることが面倒そうな式\(\ref{t_confint}\) だが，幸いにも\(t\) 分布は\(n\) が十分大きい（>50）時には正規分布と実用上問題ない程度に一致するという性質がある．正規分布の97.5%タイル点（-2.5%タイル点）は1.96だから，\(n\) が十分大きい場合は
\[
\bar{d} \pm 1.96 \cdot s_e
\]
としてスカラーデータ平均値の95%信頼区間を求めることができる．
:問題3|
スカラーデータ平均値の95%信頼区間を求めるプログラムを書いてみよう．計算のためのデータは前回同様&ref(data1.txt);を利用すること．
:問題3の解答例|
Scilabを使っている場合は&ref(calc_stats.sce);を参照．もう少しまじめなプログラム例として，&ref(calc_stats2.sci);や&ref(calc_stats3.sci);がある．
:応用問題1|
t分布を正しく計算し，データセットが小さくても正しい区間推定が行われるようなプログラムを書け．
:応用問題2|
&ref(data1.txt);はある堆積物の粒子サイズを\(\phi\) スケールで記録したものである．この堆積物の粒子の個数分布を表すヒストグラムを描け（ヒント：if文を使うことになる）．
:応用問題3|
通常，粒度分布と言ったときには粒子の個数の分布ではなく，ある粒径範囲内に含まれる粒子の総重量（総体積）の分布を意味している．画像解析で得られた粒子データの場合，面積が粒子の重量（体積）に相当するものとなる．&ref(data1.txt);の粒度分布を描き，個数分布の場合とどのような違いが見られるか観察せよ．
*角度データの平均値とベクトル集中度 [#q760e617]
角度データにはスカラーデータにはない特徴がある．それは，単純な算術平均が平均的な方向を表さないということである．例えば，0°と359°を足して2で割ると179.5°になるが，この方向はまったく平均的な方向とはいえない．
そこで，角度データの場合，それぞれの数値（角度）を方向ベクトルに変換し，すべてのデータの合成ベクトルの示す方向を平均方向とみなすことが行われてきた（図1）．
#ref(img002.png,center,zoom,320x640,図1,図1.From Borradaile (2003));
平均方向\(\bar{\theta}\) は以下の手順で計算することができる：
まず，合成ベクトルの\(x\) 成分と\(y\) 成分である\(C\) と\(S\) ，そして合成ベクトル長\(R\) を求める:
\[
\begin{eqnarray}
C = \sum{ \cos{\theta_i} },\\
S = \sum{ \sin{\theta_i} },\\
R = \sqrt{C^2 + S^2},\\
\bar{R} = R / n.
\end{eqnarray}
\]
\(\bar{R}\) はベクトル集中度（mean resultant length）と呼ばれる量で，合成ベクトル長をサンプル数で規格化したものである．\(\bar{R}\) が大きいほどデータセットが１方向にそろっていることを意味している（最小は0，最大は1）．
次に，平均方向を\(\bar{\theta}\) とすると，以下の式が成り立つ：
\[
\begin{eqnarray}
\cos{\bar{\theta}} = \frac{C}{R},\\
\sin{\bar{\theta}} = \frac{S}{R},\\
\tan{\bar{\theta}} = \frac{\sin{\bar{\theta}}}{\cos{\bar{\theta}}} = \frac{S}{C}.
\end{eqnarray}
\]
この結果，\(\bar{\theta} = \arctan{(S/C)}\) を計算すれば平均角度を求めることができそうに思える．しかし，ここで気をつけなければいけないのは，合成ベクトルがどの象限に入っているかということである．単なるアークタンジェントの計算では，\(C\) と\(S\) が共に負の場合も共に正の場合も同じ平均方向が計算されてしまう．したがって，下記のような場合わけが必要になる：
\[
\bar{\theta} = \begin{cases} \arctan{\frac{S}{C}} & (S > 0, C > 0) \\ \arctan{\frac{S}{C}}+\pi&(C < 0)\\ \arctan{\frac{S}{C}} + 2 \pi & (S < 0, C > 0)\\ \end{cases}
\]
この場合わけは，プログラミング言語によっては自動的にやってくれることもある（ExcelやJavaにはatan2，Scilabの場合はatan(y,x)）．
:問題4|
角度データの平均値とベクトル集中度（Mean Resultant Length）を求めるプログラムを書いてみよう．計算のためのデータは&ref(angledata.txt);を利用する．角度をラジアンに変換するのを忘れないこと．
:問題4の解答例|
Scilabを使っている場合は&ref(calc_cstats_avemrl.sci);を参照．
**双方向性の角度データの取り扱い [#b7462513]
堆積学では，しばしば双方向性の角度データを取り扱う．例えば，粒子の長軸方向や，グルーブキャストの方向は，双方向性のデータといえる．このようなデータを取り扱う場合,\(\theta\)と\(\theta+\pi\)が実際には同じ方向を示していることに注意しなくてはならない．
この問題を解決するには，すべての角度データをいったん2倍にして統計処理を行うと良い．最終的に，計算された平均方向を2で割れば，適切な平均方向を得ることができる．次の節で計算する信頼区間も同様である．
:問題5|
実は，先ほどの問題で扱った角度データ&ref(angledata.txt);は双方向性のものである．双方向性の角度データの平均値とベクトル集中度（Mean Resultant Length）を求めるプログラムを書いてみよう．
:問題5の解答例|
Scilabを使っている場合は&ref(calc_cstats_bd.sci);を参照．
*角度データ平均値の95%信頼区間 [#ec28e24b]
角度データセットの区間推定を行う際に問題となるのは，角度データが一般のスカラーデータセットが示すような正規分布を示さないということである．\(t\) 分布も区間推定に用いることはできない．
代わりに，比較的自然界の角度データセットに適合する理論分布として，von Mises分布と呼ばれる分布形が知られている．この分布形の性質の詳細にはここでは立ち入らない（詳しくは上記の教科書を参照）が，ここではvon Mises分布を利用した区間推定の方法について述べる．
近似解ではあるが，Fisher (1993) は比較的ベクトル集中度の高い角度データセットに対して，平均方向の標準誤差が以下の式で与えられることを示した：
\[
\frac{1}{\sqrt{ n \bar{R} \kappa}}
\]
すなわち，角度データの95%区間推定は以下のように行うことができる：
\[
\bar{\theta} \pm \arcsin{ \left( \frac{1.96}{\sqrt{ n \bar{R} \kappa}} \right) }
\]
\(\kappa\) は母集団のvon Mises分布の集中度を表すパラメーターであり，以下の方程式の解としてデータセットから推定することができる：
\[
\label{kappa_def}
I_1(\kappa)/I_0(\kappa)=\bar{R}
\]
ここで，\(I_1(x)\) と\(I_0(x)\) はそれぞれ第１次，第0次の第一種変形ベッセル関数である．
式\(\ref{kappa_def}\) は数値的に解くしかないため，\(\kappa\) を計算で求めることは非常に煩雑である．そのため，簡便な方法として以下のような表を利用することが行われている（図2）．
#ref(img003.png,center,zoom,320x640,図2,図2.From Davis (2002));
この表を用いれば，\(\bar{R}\) に対応する\(\kappa\) の推定値を簡単に読み取ることができる．この演習では，この表の数値を取り込んだテキストファイルを用意した（&ref(kappa_table.txt);）ので，プログラム制作に役立てて欲しい．
:問題6|
角度データの平均値，ベクトル集中度（Mean Resultant Length），そして平均方向の95%信頼区間を求めるプログラムを書いてみよう．計算のためのデータは&ref(angledata.txt);を利用する．
:問題6の解答例|
Scilabを使っている場合は&ref(calc_cstats.sci);を参照．
[[堆積学特別演習-プログラミングの堆積学への応用]]
