PSP:Personal Software Process – Lesson6

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

Lesson5で開発したプログラムをもとに、Lesson6用の要件を追加します。

Lesson6: 「t分布関数を0.0からxの範囲で積分した結果がpとなるようなxを求める」
● プログラム仕様スクリプト
ASGKIT PROG6を参照のこと

● 実行結果
正しい積分結果値よりxをを絞り込みながら調整しつつ、シンプソン法の最適幅、積分結果、t分布関数値を得る

OK! [W=2]
Result P: 0.0
NG! [W=2] 0.062143
NG! [W=4] 0.000017
OK! [W=8]
Result P: 0.178335
NG! [W=2] 0.116564
NG! [W=4] 0.000061
OK! [W=8]
Result P: 0.31305
NG! [W=2] 0.145412
NG! [W=4] 0.00164
OK! [W=8]
Result P: 0.396
NG! [W=2] 0.139119
NG! [W=4] 0.008551
NG! [W=8] 0.000075
OK! [W=16]
Result P: 0.441942
NG! [W=2] 0.10295
NG! [W=4] 0.021035
NG! [W=8] 0.000398
OK! [W=16]
Result P: 0.466616
NG! [W=2] 0.048413
NG! [W=4] 0.035425
NG! [W=8] 0.001493
OK! [W=16]
Result P: 0.480029
NG! [W=2] 0.014094
NG! [W=4] 0.047057
NG! [W=8] 0.003926
NG! [W=16] 0.00001
OK! [W=32]
Result P: 0.487552
NG! [W=2] 0.077696
NG! [W=4] 0.052655
NG! [W=8] 0.007991
NG! [W=16] 0.00004
OK! [W=32]
Result P: 0.491935
NG! [W=2] 0.13875
NG! [W=4] 0.050897
NG! [W=8] 0.013594
NG! [W=16] 0.000132
OK! [W=32]
Result P: 0.494589
NG! [W=2] 0.195789
NG! [W=4] 0.041945
NG! [W=8] 0.020339
NG! [W=16] 0.000344
OK! [W=32]
Result P: 0.496255
NG! [W=2] 0.13875
NG! [W=4] 0.050897
NG! [W=8] 0.013594
NG! [W=16] 0.000132
OK! [W=32]
Result P: 0.494589
NG! [W=2] 0.167808
NG! [W=4] 0.04727
NG! [W=8] 0.016859
NG! [W=16] 0.000218
OK! [W=32]
Result P: 0.495514
NG! [W=2] 0.153413
NG! [W=4] 0.049303
NG! [W=8] 0.015197
NG! [W=16] 0.000171
OK! [W=32]
Result P: 0.495078
NG! [W=2] 0.146114
NG! [W=4] 0.050156
NG! [W=8] 0.014388
NG! [W=16] 0.000151
OK! [W=32]
Result P: 0.49484
NG! [W=2] 0.149772
NG! [W=4] 0.049743
NG! [W=8] 0.01479
NG! [W=16] 0.000161
OK! [W=32]
Result P: 0.494961
NG! [W=2] 0.151592
NG! [W=4] 0.049525
NG! [W=8] 0.014993
NG! [W=16] 0.000165
OK! [W=32]
Result P: 0.49502
NG! [W=2] 0.150681
NG! [W=4] 0.049634
NG! [W=8] 0.014892
NG! [W=16] 0.000162
OK! [W=32]
Result P: 0.49499
NG! [W=2] 0.151136
NG! [W=4] 0.049578
NG! [W=8] 0.014943
NG! [W=16] 0.000164
OK! [W=32]
Result P: 0.495005
NG! [W=2] 0.150912
NG! [W=4] 0.049609
NG! [W=8] 0.014916
NG! [W=16] 0.000164
OK! [W=32]
Result P: 0.494998
NG! [W=2] 0.151025
NG! [W=4] 0.049595
NG! [W=8] 0.014929
NG! [W=16] 0.000163
OK! [W=32]
Result P: 0.495001
NG! [W=2] 0.150967
NG! [W=4] 0.049602
NG! [W=8] 0.014923
NG! [W=16] 0.000163
OK! [W=32]
Result P: 0.495
Result X: 4.60400390625

● ソースコード

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

/**
 * PSP Lesson6:
 * t分布関数を0.0からxの範囲で積分した結果がpとなるようなxを求める
 */

public class Program6 {

    private double x;
    private int dof;
    private double p;
    private double d;
    private int W;
    private int initW;
    private double E;
    private int checkF1;
    private int checkF2;
    private long adjustCount;

    public Program6() {
        initW = 1;
        W = initW;
        E = 0.00001;
        d = 0.5;
        checkF1 = 0;
        checkF2 = 0;
        adjustCount = 0;
    }

    public static void main(String[] args) {

        Program6 pg6 = new Program6();

        Double doubleArray[] = new Double[2];

        // 実行時引数からx(trial value), dof, pを取得する
        pg6.x = Double.parseDouble(args[0]);
        pg6.dof = Integer.parseInt(args[1]);
        pg6.p = Double.parseDouble(args[2]);
        boolean loopF = true;

        while (loopF) {

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

                doubleArray[i] = (double) 0;

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

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

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

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

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

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

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

                    doubleArray[i] += p6;

                }

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

                if (i == 0) {
                    pg6.W = pg6.W * 2;
                } else {
                    if ((doubleArray[i - 1] - doubleArray[i]) < pg6.E
                            && (doubleArray[i] - doubleArray[i - 1]) < pg6.E) {
                        System.out.println("OK! [W=" + pg6.W + "]");
                        System.out.println("Result P: " + doubleArray[i]);

                        if (doubleArray[i] != pg6.p) {
                            pg6.adjustD(doubleArray[i]);
                            pg6.x += pg6.d;
                            pg6.W = pg6.initW;
                            i = -1;
                        } else {
                            System.out.println("Result X: " + pg6.x);
                            loopF = false;
                        }
                    } else {
                        if (doubleArray[i - 1] > doubleArray[i]) {
                            System.out.println("NG! [W="
                                    + pg6.W
                                    + "] "
                                    + new DecimalFormat("0.######")
                                            .format(doubleArray[i - 1]
                                                    - doubleArray[i]));
                        } else {
                            System.out.println("NG! [W="
                                    + pg6.W
                                    + "] "
                                    + new DecimalFormat("0.######")
                                            .format(doubleArray[i]
                                                    - doubleArray[i - 1]));
                        }
                        // 結果保持
                        doubleArray[i - 1] = doubleArray[i];

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

            }
        }
    }

    /**
     * dを調整する
     * 
     */
    private void adjustD(Double result) {

        adjustCount += 1;
        double w = (double) (result - p);

        if (d < 0) {
            d = -d;
        }
        if (adjustCount == 1) {
            if (w > 0) {
                checkF1 = -1;
                checkF2 = -1;
            } else {
                checkF1 = 1;
                checkF2 = 1;
            }
        }

        // 正しい積分結果値(P)を通り越えるまでは、0.5ポイント単位で増減させる
        // 通過後は、より詳細に増減させて積分結果値(P)に近づくように調整する
        if (checkF1 == checkF2) {
            if (w > 0) {
                d = -0.5;
                checkF2 = -1;
            } else {
                d = 0.5;
                checkF2 = 1;
            }
        } else {
            if (w > 0) {
                d /= 2;
                d = -d;
            } else {
                d /= 2;
            }
        }

    }

    /**
     * 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();

    }

}