重回帰式の説明変数を自動選択するプログラム

このページでは「重回帰式の説明変数を自動選択するプログラム」を紹介しています。

コード

重回帰式の説明変数を自動で選択するプログラムのコードです。言語はJavaを使用しています。 統計の知識に関しては、『多変量解析法入門 (ライブラリ新数学大系) 』(永田靖)の第5章 重回帰分析を参考にしています。

計算クラス(Calculator.java)

計算を行うクラスです。

import java.util.ArrayList;
import java.util.List;

/**
 * 計算クラス
 */
public class Calculator {

	/**
	 * 説明変数を選択するかどうか?
	 * @param se1 モデル1の残差平方和
	 * @param faie1 モデル1の自由度
	 * @param se2 モデル2の残差平方和
	 * @param faie2 モデル2の自由度
	 * @return
	 */
	public boolean isSelectExplanatoryVariables(Double se1, int faie1, Double se2, int faie2) {
		Double f0 = valueForSelectExplanatoryVariables(se1, faie1, se2, faie2);
		if (f0 > 2) {
			return true;
		}
		return false;
	}

	/**
	 * 説明変数を選択するためのF0値を計算する
	 * @param se1 モデル1の残差平方和
	 * @param faie1 モデル1の自由度
	 * @param se2 モデル2の残差平方和
	 * @param faie2 モデル2の自由度
	 * @return
	 */
	public Double valueForSelectExplanatoryVariables(Double se1, int faie1, Double se2, int faie2) {
		return ((se1 - se2) / (faie1 - faie2)) / (se2 / faie2);
	}

	/**
	 * 説明変数がn個の場合の自由度を計算する
	 * @param itemXin 項目リスト(Xip)
	 * @param itemYi 項目リスト(Yi)
	 * @return 結果
	 */
	public int degreeOfFreedomN(
			List<List<Double>> itemsXip, final List<Double> itemsYi) {
		int n = itemsYi.size();
		int p = itemsXip.size();
		return n - p -1;
	}


	/**
	 * 説明変数がn個の場合の自由度調整済寄与率を計算する
	 * @param itemXin 項目リスト(Xip)
	 * @param itemYi 項目リスト(Yi)
	 * @return 結果
	 */
	public Double degreeOfFreedomAdjustedContributionRateN(
			List<List<Double>> itemsXip, final List<Double> itemsYi) {
		int n = itemsYi.size();
		int p = itemsXip.size();
		int faiT = n - 1;
		int faie = n - p -1;
		return 1 - ((residualSumOfSquaresN(itemsXip, itemsYi) / faie) / (sumOfSquares(itemsYi) / faiT));
	}

	/**
	 * 説明変数がn個の場合の残差平方和を計算する
	 * @param itemXin 項目リスト(Xip)
	 * @param itemYi 項目リスト(Yi)
	 * @return 結果
	 */
	public Double residualSumOfSquaresN(
			List<List<Double>> itemsXip, final List<Double> itemsYi) {
		return sumOfSquares(itemsYi) - multipleCorrelationCoefficient(itemsXip, itemsYi);
	}

	/**
	 * 説明変数がn個の場合の回帰による平方和を計算する
	 * @param itemXin 項目リスト(Xip)
	 * @param itemYi 項目リスト(Yi)
	 * @return 結果
	 */
	public Double multipleCorrelationCoefficient(
			List<List<Double>> itemsXip, final List<Double> itemsYi) {
		List<Double> slist = new ArrayList<>();

		for (List<Double> item : itemsXip) {
			slist.add(deviationSumOfProduct(item, itemsYi));
		}
		List<Double> batas = partialRegressionCoefficientsForExpVarN(itemsXip, itemsYi);

		double result = 0;

		for (int i = 0; i < batas.size(); i++) {
			result += (batas.get(i) * slist.get(i));
		}
		return result;
	}

	/**
	 * 説明変数がn個の場合の切片を計算する
	 * @param itemsXip 項目リスト(Xip)
	 * @param itemsYi 項目リスト(Yi)
	 * @return 結果
	 */
	public Double sectionForExpVarN(
			final List<List<Double>> itemsXip, final List<Double> itemsYi) {
		List<Double> xbars = new ArrayList<>();

		for (List<Double> item : itemsXip) {
			xbars.add(average(item));
		}
		Double ybar = average(itemsYi);
		List<Double> batas = partialRegressionCoefficientsForExpVarN(itemsXip, itemsYi);

		double result = ybar;

		for (int i = 0; i < batas.size(); i++) {
			result -= (batas.get(i) * xbars.get(i));
		}
		return result;
	}

	/**
	 * 説明変数がn個の場合の偏回帰係数を計算する
	 * @param itemsXip 項目リスト(Xip)
	 * @param itemsYi 項目リスト(Yi)
	 * @return 結果
	 */
	public List<Double> partialRegressionCoefficientsForExpVarN(
			final List<List<Double>> itemsXip, final List<Double> itemsYi) {
		List<List<Double>> coefficientMatrix = new ArrayList<>();
		int n = itemsXip.size();

		for (int i = 0; i < n; i++) {
			List<Double> row = new ArrayList<>();

			for (int j = 0; j < n; j++) {

				if (i == j) {
					row.add(sumOfSquares(itemsXip.get(i)));

				} else {
					row.add(deviationSumOfProduct(itemsXip.get(i), itemsXip.get(j)));
				}
			}
			row.add(deviationSumOfProduct(itemsXip.get(i), itemsYi));
			coefficientMatrix.add(row);
		}
		return systemOfEquations(coefficientMatrix);
	}

	/**
	 * 偏差積和を計算する
	 * @param itemsXi 項目リスト(Xi)
	 * @param itemsYi 項目リスト(Yi)
	 * @return 結果
	 */
	public Double deviationSumOfProduct(final List<Double> itemsXi, final List<Double> itemsYi) {
		List<Double> itemsXiYi = new ArrayList<>();
		int n = itemsXi.size();

		for (int i = 0; i < n; i++) {
			itemsXiYi.add(itemsXi.get(i) * itemsYi.get(i));
		}
		Double xiyiSum = sum(itemsXiYi);
		Double xiSum = sum(itemsXi);
		Double yiSum = sum(itemsYi);
		return xiyiSum - ((xiSum * yiSum) / n);
	}

	/**
	 * 平方和を計算する
	 * @param items 項目リスト
	 * @return 結果
	 */
	public Double sumOfSquares(final List<Double> items) {
		Double xbar = average(items);
		List<Double> squares = new ArrayList<>();

		for (Double item : items) {
			Double sqare = (item - xbar) * (item - xbar);
			squares.add(sqare);
		}
		return sum(squares);
	}

	/**
	 * 平均値を計算する
	 * @param items 項目リスト
	 * @return 結果
	 */
	public Double average(final List<Double> items) {
		return sum(items) / items.size();
	}

	/**
	 * 総和を計算する
	 * @param items 項目リスト
	 * @return 結果
	 */
	public Double sum(final List<Double> items) {
		Double result = 0.0;

		for (Double item : items) {
			result += item;
		}
		return result;
	}

	/**
	 * 連立方程式を計算する
	 * @param coefficientMatrix m行n列の係数行列
	 * @return 結果
	 */
	public List<Double> systemOfEquations(final List<List<Double>> coefficientMatrix) {
		int n = coefficientMatrix.size();
		int m = coefficientMatrix.get(0).size();

		double[][] matrix = new double[n][];

		for (int i = 0; i < n; i++) {
			double[] row = new double[m];

			for (int j = 0; j < m; j++) {
				row[j] = coefficientMatrix.get(i).get(j);
			}
			matrix[i] = row;
		}
		return systemOfEquations(matrix);
	}

	/**
	 * 連立方程式を計算する
	 * @param coefficientMatrix m行n列の係数行列
	 * @return 結果
	 */
	public List<Double> systemOfEquations(final double[][] coefficientMatrix) {
		List<Double> result = new ArrayList<>();
		double[][] matrix = coefficientMatrix;
		int n = matrix.length;
		int m = matrix[0].length;

		for (int k = 0; k < n; k++) {
			double p = matrix[k][k];

			for (int j = k; j < m; j++) {
				matrix[k][j] = matrix[k][j] / p;
			}

			for (int i = 0; i < n; i++) {
				if (i != k) {
					double d = matrix[i][k];

					for (int j = k; j < m; j++) {
						matrix[i][j] = matrix[i][j] - d * matrix[k][j];
					}
				}
			}
		}

		for (int i = 0; i < n; i++) {
			result.add(matrix[i][n]);
		}

		return result;
	}
}

動作確認用クラス

動作確認用のクラスです。

import java.util.ArrayList;
import java.util.List;

/**
 * サンプルのテスト用クラス
 */
public class SampleTest {

	public static void main(String[] args) {
		testExp();
	}

	static void testExp() {
		double[][] xip = {
			{51, 38, 57, 51, 53, 77, 63, 69, 72, 73},
			{16, 4, 16, 11, 4, 22, 5, 5, 2, 1},
			{1, 2, 1, 2, 1, 3, 4, 2, 2, 1}
		};

		double[] yi = {
			3.0, 3.2, 3.3, 3.9, 4.4, 4.5, 4.5, 5.4, 5.4, 6.0
		};
		selectExplanatoryVariables(xip, yi);
	}

	static void selectExplanatoryVariables(double[][] xip, double[] yi) {
		List<Double> itemsYi = prepareTestData(yi);
		List<List<Double>> itemsXipM1 = new ArrayList<>();

		Calculator calc = new Calculator();
		Double se1 = calc.sumOfSquares(itemsYi);
		int faie1 = itemsYi.size() - 1;

		for (int i = 0; i < xip.length; i++) {
			List<List<Double>> itemsXipM2 = new ArrayList<>();
			for (List<Double> item : itemsXipM1) {
				itemsXipM2.add(item);
			}
			List<Double> additional = prepareTestData(xip[i]);
			itemsXipM2.add(additional);
			Double se2 = calc.residualSumOfSquaresN(itemsXipM2, itemsYi);
			int faie2 = calc.degreeOfFreedomN(itemsXipM2, itemsYi);

			System.out.println("F0="
				 + calc.valueForSelectExplanatoryVariables(se1, faie1, se2, faie2));

			if (calc.isSelectExplanatoryVariables(se1, faie1, se2, faie2)) {
				itemsXipM1.add(additional);
			}
			se1 = calc.residualSumOfSquaresN(itemsXipM1, itemsYi);
			faie1 = calc.degreeOfFreedomN(itemsXipM1, itemsYi);
		}
		printResult(itemsXipM1, itemsYi);
	}

	static void printResult(List<List<Double>> itemsXip, List<Double> itemsYi) {
		Calculator calc = new Calculator();
		Double beta0 = calc.sectionForExpVarN(itemsXip, itemsYi);
		List<Double> batas = 
			calc.partialRegressionCoefficientsForExpVarN(itemsXip, itemsYi);

		int n = batas.size();
		StringBuilder msg1 = new StringBuilder();
		StringBuilder msg2 = new StringBuilder();
		String fmt1 = "切片(β0)\t\t:%f\n";
		String fmt2 = "偏回帰係数(β%d)\t:%f\n";
		String fmt3 = "重回帰式(p=%d)\t\t:y = %s\n";
		msg1.append(String.format(fmt1, beta0));
		msg2.append(beta0);

		for (int i = 0; i < n; i++) {
			msg1.append(String.format(fmt2, i + 1, batas.get(i)));
			msg2.append(getMark(batas.get(i)) + batas.get(i) + "x" + i + 1);
		}
		System.out.println(
			msg1.toString() + String.format(fmt3, n, msg2.toString()));
		System.out.println("重相関係数:"
			 + calc.multipleCorrelationCoefficient(itemsXip, itemsYi));
		System.out.println("自由度調整済寄与率:"
			 + calc.degreeOfFreedomAdjustedContributionRateN(itemsXip, itemsYi));
	}

	static List<Double> prepareTestData(double[] sample) {
		List<Double> items = new ArrayList<>();

		for (double data : sample) {
			items.add((double) data);
		}
		return items;
	}

	static String getMark(Double data) {
		String mark = " + ";
		if (data < 0) {
			mark = " ";
		}
		return mark;
	}
}

実行結果

実行すると以下のように表示されます。F0の値が2以上を基準に変数が選択されます。

F0=13.439322860158985
F0=43.59889285283754
F0=0.4205569318080352
切片(β0)		:1.020130
偏回帰係数(β1)	:0.066805
偏回帰係数(β2)	:-0.080830
重回帰式(p=2)		:y = 1.020129547203318 + 0.06680476622865743x01 -0.0808299334202589x11

重相関係数:8.937513015337364
自由度調整済寄与率:0.9336286856966449

戻る

Copyright (C) 2016 計算プログラムのコード. All Rights Reserved. Loarding…