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();
}
}
に沿いながら、8本ほどのプログラムを開発し、開発プロセスの見直しをしました。