行列の統計処理の例
この節では、JSLによる行列を使った統計処理の例を紹介します。
回帰の例
回帰分析の計算手法について、JMPに組み込まれている機能を使わず、ユーザ自身が独自の方法を実施したい場合を考えてみましょう。そのような場合、簡潔な行列表現を使って、ほんの数行で書くことができます。
Y = [98, 112.5, 84, 102.5, 102.5, 50.5, 90, 77, 112, 150, 128, 133, 85, 112];
X = [65.3, 69, 56.5, 62.8, 63.5, 51.3, 64.3, 56.3, 66.5, 72, 64.8, 67, 57.5, 66.5];
X = J( N Row( X ), 1 ) || X; // 切片の列に1を入れる
beta = Inv( X` * X ) * X` * Y; // 最小2乗推定
resid = Y - X * beta; // 残差、Y - 推定値
sse = resid` * resid; // 誤差平方和
Show( beta, sse );
これは、データテーブルからデータを取得し、関連する統計量を計算してレポートウィンドウに表示するスクリプトへと拡張できます。
dt = Open( "$SAMPLE_DATA/Big Class.jmp" );
// データを行列に移す
x = (Column( "年齢" ) << Get Values ) || (Column( "身長(インチ)" ) << Get Values);
x = J(N Row( x ), 1, 1) || x;
y = Column( "体重(ポンド)" ) << Get Values;
// 回帰分析の計算
xpxi = Inv( x` * x );
beta = xpxi * x` * y; // パラメータ推定値
resid = y - x * beta; // 残差
sse = resid` * resid; // 誤差平方和
dfe = N Row( x ) - N Col( x ); // 自由度
mse = sse / dfe; // 誤差の平均平方、誤差分散推定値
// 推定値の追加計算
stdb = Sqrt( Vec Diag( xpxi ) * mse ); // 推定値の標準誤差
alpha = .05;
qt = Students t Quantile( 1 - alpha / 2, dfe );
betau95 = beta + qt * stdb; // 上側95％信頼限界
betal95 = beta - qt * stdb; // 下側95％信頼限界
tratio = beta :/ stdb; // Studentのt値
probt = ( 1 - TDistribution( Abs( tratio ), dfe ) ) * 2; // p値
// 結果の表示
New Window( "ビッグ クラスの回帰分析",
Table Box(
String Col Box( "項",{"切片", "年齢", "身長(インチ)"} ),
Number Col Box( "推定値", beta ),
Number Col Box( "標準誤差", stdb ),
Number Col Box( "t値", tratio ),
Number Col Box( "p値", probt ),
Number Col Box( "下側95%", betal95 ),
Number Col Box( "上側95%", betau95 ) ) );
分散分析（ANOVA）の例
行列演算を用いて一元配置分散分析も実行できます。この例では、3水準の要因（薬品の投与量で低用量、中用量、高用量）と1つの応答変数から成る問題を扱います。つまりこの例は、次のような一般線形モデルで表されます。
Y = a + bX + e
ここで
• Yは応答変数のベクトル
• aは切片項
• bは係数から成るベクトル
• Xは要因の計画行列
• eは誤差項
factor = [1, 2, 3, 1, 2, 3, 1, 2, 3];
y = [1, 2, 3, 4, 3, 2, 5, 4, 3];
まず、要因の計画行列を作成する必要があります。
Design Nom( factor );
[1 0, 0 1, -1 -1, 1 0, 0 1, -1 -1, 1 0, 0 1, -1 -1]
次に、切片項として、計画行列に1だけから成る列を加えます。これを行うには、単にJとDesign Nom()を結合するだけです。
x = J( 9, 1, 1) || Design Nom( factor );
[1 1 0,
1 0 1,
1 -1 -1,
1 1 0,
1 0 1,
1 -1 -1,
1 1 0,
1 0 1,
1 -1 -1]
ここで、一般的な方程式を解くために、次のような行列Mを作成する必要があります。
行列Mは、次のような結合をすれば1ステップで作成できます。
M = ( x` * x || x` * y )
|/
( y` * x || y` * y );
[9 0 0 27,
0 6 3 2,
0 3 6 1,
27 2 1 93]
モデル全体の推定結果は、行列Mにおいて、X′Xの部分をすべて掃き出せば得られます。また、最初の列だけ掃き出せば切片のみのモデルの推定結果が得られます。
FullFit = Sweep( M, [1, 2, 3] ); // フルモデルのあてはめ
InterceptOnly = Sweep( M, [1] ); // 切片のみのモデル
分散分析の一部の結果は、これら2つのモデルを比較して計算されます。ここでは、フルモデル（切片と傾きを含んだモデル）の結果を見てみましょう。フルモデルに関して掃き出された行列は、次のような行列になります。
[0.111111111111111 0 0 3,
0 0.222222222222222 -0.11111111111111 0.333333333333333,
0 -0.11111111111111 0.222222222222222 0,
-3 -0.33333333333333 0 11.33333333333333]
行列右上部のモデル係数を見てください。左下部の係数も、符号が正負逆になっている以外は同じ、3、0.333、0です。結果は次のように解釈できます。
• 切片項の係数は3です。
• 要因の1番目の水準の係数は0.333です。
• 2番目の水準の係数は0です。
• Design Nom()を使っているので、3番目の水準の係数は-0.333です。
• 行列右下部には誤差平方和 11.333が入ります。
このプログラム例では、特定の数値を指定していますが、それらを適切な引数に置き換えることにより、より汎用性があるスクリプトに変更できます。上の結果は、モデルのあてはめプラットフォームの結果と一致しています。
図7.1 モデルのあてはめのANOVAレポート
Figure 7.1のレポートは次のようにして作成します。
1. データテーブルで説明している方法で、データテーブルを作成します。
dt = New Table( "Data" );
dt << New Column( "y", Set Values( [1, 2, 3, 4, 3, 2, 5, 4, 3] ) );
dt << New Column( "factor", "Nominal", Values( [1, 2, 3, 1, 2, 3, 1, 2, 3] ) );
2. 分析を起動し、そのスクリプトをウィンドウに保存するで説明している方法で、モデルのあてはめを実行します。
obj = Fit Model(Y( :y ),
Effects( :factor ),
Personality( "Standard Least Squares" ),Run(
y << {Plot Actual by Predicted( 0 ), Plot Residual by Predicted( 0 ), Plot Effect Leverage( 0 )})
);
3. ディスプレイボックスオブジェクトの参照の取得で説明している方法で、ディスプレイを操作します。
ranova = obj << Report;ranova[Outline Box( 6 )] << Close( 0 );
ranova[Outline Box( 7 )] << Close( 1 );
ranova[Outline Box( 9 )] << Close( 1 );
ranova[Outline Box( 5 ), Number Col Box( 2 )] << Select;
ranova[Outline Box( 6 ), Number Col Box( 1 )] << Select;
この結果がFigure 7.1に示されているものになります。