SASユーザー会投稿論文2006

SAS を用いた数値計算および天体データの統計処理
作花一志
南野公彦
京都情報大学院大学
Educational Programs of Numerical Computation
and Statistical Study of Some Celestial Bodies
by using
Kazushi Sakka
and
SAS
Kimihiko Nouno
The Kyoto College of Graduate Studies for Informatics
数値計算・統計解析・多変量解析の数学的基礎を学ぶための SAS プロシ―ジャを用
いた教材と SAS の天文学への応用例を紹介する。このためには SAS/iml のプロシージ
ャを用いると効果的であることがわかった。また自然科学データの統計処理の例として
1100 個のカイパーベルト天体についての結果を報告する。
なおこれは 7 月 27 日に SAS ユーザー会学術総会で発表した論文を一部改訂したもの
である
グラフ描画
PHP
非線形回帰分析
クラスター分析
カイパーベルト天体
１．SAS による数値計算
1.1
SAS による数値計算
数値計算はコンピュータの基本的かつ重要な利用シーンの一つである。しかし，近年の
初・中等教育機関での数学の履修時間の減少の影響か，数値計算法を（本当はそのバック
ボーンである数学そのものを）学生に理解させるのは困難になりつつある。
一方で，プログラミング言語（C，Fortran，Java など）では既に膨大な数の数値計算用
ライブラリが開発されているし，数値計算用パッケージ（MATLAB，Octave，Scilab など），
データ解析用パッケージ（SAS，SPSS，S/S-PLUS，R など）も非常に広範囲の数値計算
法をカバーしている。ブラックボックスであることを認めてしまえば，プログラムを一か
ら作るよりもこれらのパッケージを活用したほうが，開発効率，計算速度，バグの少なさ
など，どれをとっても優れている。
そこで，筆者らは数値計算の入門的なアルゴリズムの仕組みをグラフに表示し，学生に
対してビジュアル的に訴えることで，数学に不慣れな学生にも数値計算がどのようなもの
か，パッケージが内部で何を行っているのかを体感してもらうことだけでも達成したいと
考えて教材の作成を行っている。その際，行列演算を容易に表現でき，数学関数・サブル
ーチンも多数用意されている SAS/IML と，強力なグラフ描画力を持つ SAS/GRAPH の機
能を使用した。
1.2
SAS ソースコードを作成する PHP プログラミング
SAS プログラミング作成に当たって，エディタの上でコードを書くのではなく，ウェブ
から必要なパラメータを入力して作成する方法を開発中である。現在，次の 6 個について
は作成済みであり，筆者の Web ページ（http://web1.kcg.edu/~sakka/PHP/calc/）から利用
できる。開発には PHP を使用したが，閲覧は IE で可能である。
a) 整方程式の解法
b) 非線型方程式の解法（ニュートンラフソン法）
c) 数値積分（シンプソン則）
d) 極座標表示グラフ描写
e) 立体図形グラフ描写
f) 行列演算（行列式
固有値など）
図 1 は立体図形グラフを描写する場合で，関数 z＝f(x,y) を選ぶかまたは記すことでサー
バに SAS コードが生成され，テキストボックスに表示される。それをダウンロードしてロ
ーカルフォルダーに保存して，SAS 9.1 から起動・実行することによって，図 2，図 3 のよ
うな結果が得られる。ここで使用するプロシージャは proc g3d である。
図 1 立体図形グラフ作成パラメータ入力画面
図 2 アウトプット画面
z＝x2―y2
図 3 アウトプット画面
z＝x2＋y2
また行列演算パラメータ入力画面において 9 個の数値を入力すると，
図 4 のように行列式，
固有値，固有ベクトル，逆行列が得られる。
図 4 行列演算の結果
今後は次のテーマに取り組んでいく予定である。
a) 微分方程式の解法
b) 線形計画法
c) 統計解析
集計
d) 多変量解析
回帰・主成分・クラスター分析
相関
確率分布など
この開発は，2006 年度の京都情報大学院大学（http://www.kcg.edu）の修了プロジェク
トのテーマにもなっており，その結果については，次回で報告を行いたい。
２．カイパーベルト天体の統計的性質
2.1
KBO とは
カイパーベルト天体（KBO）とは冥王星の彼方にある小天体で，1992 年に発見されて以
来，現在では約 1200 個が確認されている。サイズは地球の月よりも小さいが，地球や他の
惑星のように太陽の周りを公転している天体である。今年 8 月に惑星から除外された冥王
星もグループに含まれる。その運動は太陽からの引力によるので 3 元連立 2 回微分方程式
で表され，6 個の積分定数をもち，ある平面内で楕円軌道を描く。太陽系天体の運動を記
述する際には日心黄道座標が使用される。それは図 5 に示すように中心は太陽で，春分
点方向をｘ軸、地球の軌道面（黄道面）をｘｙ面とする座標系であり，ここでａ（平均
距離）ｅ（離心率）Ω（昇交点黄経）ω（近日点引数）i（軌道傾斜）は軌道要素と呼ば
れる各惑星固有の定数である。
図 5 惑星の軌道要素
http://cfa-www.harvard.edu/iau/lists/MPLists.html に載っている 1100 個の KBO に冥王星
を含めてそれらの軌道要素について、SAS/STAT で統計的調査を行った。絶対等級 H とは
天体の明るさを表すもので値が小さいほど明るい。
2.2 軌道要素の分布
データセットよりソリューション－分析－分布
の操作より次の諸量が得られる。
表 1 KBO の軌道要素の統計諸量
平均
平均距離 a
中央値
標準偏差
最大
最小
45.1816
43.7155
24.3827
534
5.206
離心率 e
0.1509
0.102
0.1595
0.96
0
軌道傾斜角 i
8.9165
5.200
9.0712
78.0
0.20
絶対等級 H
7.4227
7.300
1.4912
14.3
－1.2
絶対等級 H の分布はほぼ正規分布であるが，もっと遠くのもっと暗い KBO が見つかれば
右上がりの指数分布になるかもしれない。注目すべきは 3 個だけ異常に明るいものがある
ことだ。最も明るい天体はやはりエリス（昨年発見され第十惑星として話題になった
2003UB313）で次いで冥王星であった。平均距離は 2 個ほど例外的に大きいものが見つか
っているが，ほとんどは平均値あたりに分布する。
図 6 絶対等級の分布
図 7 平均距離の分布
一方，離心率 e は 0 から 1 までの，軌道傾斜 i は 0 から 180 までの数であるが，小さな値
に偏っている。
図 8 離心率の分布
2.3 軌道要素の散布図
図 9 はe とaの散布図であり，各点はシャープなエッジの下に現れ， eが小さくaが大き
いものは見られない。これは明らかに選択効果である。近日点距離ａ(１―ｅ)がある値以下
のものしか観測できないのだから
ａ(１―ｅ)＜a0 すなわち ａ＜a0／(１―ｅ)
これが包絡線を表し，a0 は約 40 と読み取れる。最もeが大きいKBOはまた最もaが大きい
ということは自明であろう。
図 9 散布図
離心率と平均距離
H と a， H と i，H と e， e と i， a と i の間に顕著な傾向は見られなかった。
2.4 非線形回帰分析
軌道傾斜角 i の分布は図 10 のようになり，
図 10 軌道傾斜角の分布
これを
y = (ax+b)exp(－x/c)
で近似し３つの係数を回帰分析より求めよう。残差平方和を最小にするため，ｙをａ，ｂ，
ｃで偏微分して 0 とおくとａ，ｂ，ｃに関する連立方程式が得られるが，非線形であるの
で，初期値を与え繰り返し演算を行う。
∂y/∂a＝x･exp(－x/c)
∂y/∂b＝exp(－x/c)
∂y/∂c＝(ax+b )exp(－x/c)( x/c2)
非線形回帰の場合は reg プロシージャではなく nlin プロシージャを使う。
係数，予測値・上下 95%近似信頼限界・残差を求め，結果はデータセット d2 に出力する。
proc nlin data=d1;
/* パラメータの初期値 */
parms a=100 b=20 c=10;
model y=(a*x+b)*exp(－x/c);
/* 回帰式 */
output out=d2 p=pred u95m=clu95 l95m=cll95; /* データセット書き出し */
der.a= x*･exp(－x/c);
/* ∂y/∂a */
der.b= exp(－x/c);
/* ∂y/∂b
der.c=(a*x+b )*exp(－x/c)*( x/c/c);
/* ∂y/∂c */
run;
この結果 14 回の繰り返しで次の値が得られた。
表 2 非線形回帰の収束の様子
a=165.3
b=17.76
c=2.15
*/
このデータセットにてグラフを描く。
proc gplot data=d2;
plot (pred clu95 cll95 y)*x/overlay;
symbol1 v==
c=blue i=spline l=1;
/* pred を実線で */
symbol2 v=none c=red i=spline l=33;
/* clu を破線で */
symbol3 v=none c=red i=spline l=33;
/* cll を破線で */
symbol4 v=square c=magenta i= ;
/* ｙを線なしで */
run;
図 11 回帰分析の結果のグラフ
グラフで▢印は実測値，実線は予測値，点線は上下 95%近似信頼限界を表す。しかしあま
り良い近似ではなく，特に 10＜ｘ1＜30 では実測値と予測値がずれている。
この回帰分析では[1]に載っている方法を使用した。
2.5
クラスター分析
1100 個の天体には重複して登録されているものがあるかもしれない。それを調べるため
に a，e，ｉ，ω，Ω，についてクラスター分析を行った。すなわちデータセットの 5 変数の
作る 5 次元空間における距離を求め，小さい順に列挙していく。この方法で楕円のサイズ・
扁平度・傾き・方向において非常に似通った 6 個の天体が見つかった。2002TL301，
2002TB301，2002TJ301，2002TD301，2002TF301，2002TG301 の 6 個の天体は同一の
ものが重複して登録されているのではなく，現在同じ方向に見えているだけのようだが，
天文学考察は割愛する。
表 3 クラスター分析の結果
図 12 は作花・中西[2]の惑星軌道プログラム（VisualBasic）を使い，2006 年 7 月 26 日に
おける地球・天王星・海王星・冥王星とこの６つの KBO について黄道面に投影された軌道
を描いたものである。その形は a，e，pi（＝ω＋Ω）で決まり，a は楕円の大きさを e は偏
平率を pi は楕円の長軸の方向を表す。6 個の KBO はほぼ同じ方向に集合している。
図 12 KBO の軌道図
KBO の軌道要素について宇宙航空研究開発機構の吉川真氏より貴重なアドバイスを頂い
たこと厚くお礼申し上げるしだいである。
参考文献
[1] 竹内 啓 監修
[2]
「SAS で学ぶ統計的データ解析 1-7」
作花一志・中西久崇「天文学入門」オーム社
2001
東京大学出版会 1996