PSP:Personal Software Process – Lesson5

Facebook にシェア
このエントリーをはてなブックマークに追加
[`livedoor` not found]
Delicious にシェア

Lesson5: 「シンプソン法を使用して、t分布関数を数値的に分析する」
● プログラム仕様スクリプト
ASGKIT PROG5を参照のこと

● 実行結果
シンプソン法の幅をより詳細に設定していき、差分が0.00001以下であれば適切な幅とみなし終了する。

NG! [W=2] 0.094385
NG! [W=4] 0.02666
NG! [W=8] 0.000162
OK! [W=16]
Result: 0.494999

● ソースコード

import java.math.BigDecimal;
import java.text.DecimalFormat;

/**
 * PSP Lesson5:
 * シンプソン法を使用して、t分布関数を数値的に分析する
 */

public class Program5 {

    private double x;
    private int dof;
    private int W;
    private double E;

    public static void main(String[] args) {

        Program5 pg5 = new Program5();

        Double doubleArray[] = new Double[2];

        // 実行時引数からx, dof, Wを取得する
        pg5.x = Double.parseDouble(args[0]);
        pg5.dof = Integer.parseInt(args[1]);
        pg5.W = Integer.parseInt(args[2]);
        pg5.E = 0.00001;

        for (int i = 0; i < 2; i++) {

            doubleArray[i] = (double) 0;

            // wを計算する
            double w = pg5.x / pg5.W;

            for (int j = 0; j <= pg5.W; j++) {

                // p2を計算する
                double p2 = 1 + (Math.pow((double) w * j, 2) / pg5.dof);
                p2 = pg5.changeScale(p2, 6);

                // p3を計算する
                double p3 = Math.pow(p2, -(double) (pg5.dof + 1) / 2);
                p3 = pg5.changeScale(p3, 6);

                // p4を計算する
                double p4 = pg5.calcP4(p3, pg5.dof);
                p4 = pg5.changeScale(p4, 6);

                // p5を計算する
                double p5 = p3 * p4;
                p5 = pg5.changeScale(p5, 6);

                // p6を計算する
                double p6 = pg5.calcP6(w, p5, j);

                doubleArray[i] += p6;

            }

            doubleArray[i] = pg5.changeScale(doubleArray[i], 6);

            if (i == 0) {
                pg5.W = pg5.W * 2;
            } else {
                if ((doubleArray[i-1] - doubleArray[i]) < pg5.E
                        && (doubleArray[i] - doubleArray[i-1]) < pg5.E) {
                    System.out.println("OK! [W=" + pg5.W + "]");
                    System.out.println("Result: " + doubleArray[i]);
                } else {
                    if (doubleArray[i-1] > doubleArray[i]) {
                        System.out
                                .println("NG! [W="
                                        + pg5.W
                                        + "] "
                                        + new DecimalFormat("0.######")
                                                .format(doubleArray[i-1]
                                                        - doubleArray[i]));
                    } else {
                        System.out
                                .println("NG! [W="
                                        + pg5.W
                                        + "] "
                                        + new DecimalFormat("0.######")
                                                .format(doubleArray[i]
                                                        - doubleArray[i-1]));
                    }
                    // 結果保持
                    doubleArray[i - 1] = doubleArray[i];

                    // シンプソン法適用時の幅拡張
                    pg5.W = pg5.W * 2;
                    i = 0;
                }
            }
        }
    }

    /**
     * p4を計算する
     * 
     */
    private Double calcP4(Double p3, int dof) {

        double base = fact(((double) (this.dof + 1) / 2) - 1);
        double frac = fact((double) this.dof / 2 - 1.0);
        return base / (Math.pow(this.dof * Math.PI, 0.5) * frac);

    }

    /**
     * p6を計算する
     * 
     */
    private Double calcP6(Double w, Double p5, int i) {

        int multi = 0;

        if (i == 0 || i == this.W) {
            multi = 1;
            // 偶数の場合はmulti = 2
        } else if (i % 2 == 0) {
            multi = 2;
            // 奇数の場合はmulti = 4
        } else {
            multi = 4;
        }

        return (double) (w / 3) * multi * p5;

    }

    /**
     * 階乗計算
     * 
     */
    private double fact(double i) {

        if (i == 1) {
            return 1;
        } else {
            if (i == 0) {
                return i;
            } else if (i < 1) {
                return i * Math.sqrt(Math.PI);
            } else {
                // 再帰呼び出しする
                return i * fact(i - 1);
            }
        }

    }

    /**
     * 四捨五入を行う
     * 
     */
    private double changeScale(double val, int i) {

        BigDecimal bd = new BigDecimal(val);
        return bd.setScale(i, BigDecimal.ROUND_HALF_UP).doubleValue();

    }

}