diff --git a/pom.xml b/pom.xml
index 3c91035..df06e2e 100644
--- a/pom.xml
+++ b/pom.xml
@@ -56,9 +56,14 @@
jlfgr
1_0
+
+
+ org.apache.commons
+ commons-math3
+ 3.6.1
+
-
diff --git a/src/main/java/info/istlab/Zemi01/Ttest.java b/src/main/java/info/istlab/Zemi01/Ttest.java
index ac2fca1..682a4b4 100644
--- a/src/main/java/info/istlab/Zemi01/Ttest.java
+++ b/src/main/java/info/istlab/Zemi01/Ttest.java
@@ -1,22 +1,29 @@
package info.istlab.Zemi01;
+import org.apache.commons.math3.distribution.TDistribution;
+
public class Ttest {
public static void main(String[] args) {
double[] d1 = { 321, 269, 341, 318, 297, 308, 343, 337, 329 };
double[] d2 = { 303, 253, 283, 333, 307, 310, 296, 265, 271 };
double t1, t2;
System.out.println(t1 = paired_t(d1, d2));
- double p1 = tdist(t1, d1.length-1, 2);
- System.out.println(p1);
-
+ double p1 = tdist(t1, d1.length - 1, 2);
+ System.out.println("p1 = "+p1);
System.out.println("----");
- System.out.println(t2 = welch_t(d1, d2));
+ System.out.println("t2 = "+ (t2 = welch_t(d1, d2)));
double df2 = welch_degree_of_freedom(d1, d2);
- System.out.println("df2 = "+df2);
- double p2 = tdist(t2,df2, 1);
- System.out.println(p2);
+ System.out.println("df2 = " + df2);
+ double p2 = tdist(t2, df2, 2);
+ System.out.println("p2 = "+p2);
+
+ // System.out.println("----");
+ // TDistribution dist = new TDistribution(2);
+ // for(float x = -5; x < 5; x+=0.1){
+ // System.out.println(x+" "+dist.cumulativeProbability(x));
+ // }
}
/**
@@ -101,25 +108,45 @@
double diff_avg = avg(ary1) - avg(ary2);
return diff_avg / Math.sqrt(stddev1 + stddev2);
}
- public static double welch_degree_of_freedom(double[] ary1, double[] ary2){
+
+ public static double welch_degree_of_freedom(double[] ary1, double[] ary2) {
// 不偏分散/データ数
double stddev1 = var(ary1) / ary1.length;
double stddev2 = var(ary2) / ary2.length;
double a1 = Math.pow(stddev1 + stddev2, 2);
- double b1 = Math.pow(stddev1,2)/(ary1.length-1) + Math.pow(stddev2,2)/(ary2.length-1);
- return a1/b1;
+ double b1 = Math.pow(stddev1, 2) / (ary1.length - 1) + Math.pow(stddev2, 2) / (ary2.length - 1);
+ return a1 / b1;
}
/**
- * TDIST
+ * TDIST by commons-math3 TDistribution
+ *
* @param x
* @param df
* @param tails
* @return
*/
- public static double tdist(double x, double deg_free, int tails){
- return (1.0d - studentsCdf(x, deg_free)) * tails;
+ public static double tdist(double x, double deg_free, int tails) {
+ TDistribution tdistobj = new TDistribution(deg_free);
+ double cdf = tdistobj.cumulativeProbability(x);
+ double ret = (1-cdf) * tails;
+ // //test
+ // double test = tdist_old(x, deg_free, tails);
+ // if (test != ret) System.out.println("DIFF ret = "+ret+" test "+test );
+ return ret;
}
+
+ /**
+ *
+ * @param x
+ * @param deg_free
+ * @param tails
+ * @return
+ */
+ public static double tdist_old(double x, double deg_free, int tails) {
+ return (1-studentsCdf(x, deg_free)) * tails;
+ }
+
/**
*
* http://www.java2s.com/example/java-api/java/lang/math/exp-1-32.html
@@ -163,7 +190,7 @@
* Log Gamma Function
*
* @param Z
- * @return
+ * @return
*/
public static double logGamma(double Z) {
double S = 1.0 + 76.18009173 / Z - 86.50532033 / (Z + 1.0) + 24.01409822 / (Z + 2.0)
@@ -172,13 +199,14 @@
return LG;
}
+
/**
* Internal function used by StudentCdf
*
* @param x
* @param A
* @param B
- * @return
+ * @return
*/
protected static double betinc(double x, double A, double B) {
double A0 = 0.0;
@@ -187,7 +215,7 @@
double B1 = 1.0;
double M9 = 0.0;
double A2 = 0.0;
- while (Math.abs((A1 - A2) / A1) > 0.00000000000000001) {
+ while (Math.abs((A1 - A2) / A1) > 0.00000001) {
A2 = A1;
double C9 = -(A + M9) * (A + B + M9) * x / (A + 2.0 * M9) / (A + 2.0 * M9 + 1.0);
A0 = A1 + C9 * A0;
@@ -204,4 +232,78 @@
return A1 / A;
}
+ /**
+ * うまくいかないので使わない
+ * @param x
+ * @param v
+ * @return
+ */
+ public static double studentsCdf2(double x, double v) {
+ double xL = x;
+ double xR = 100000;
+ int n = 10000;
+ double goukei = 0.0;
+ double a, b, xc, tmp;
+ double dx = (xR - xL) / n;
+ double B = simpson_betafunc(0.5, v / 2);
+ // double B = Beta.regularizedBeta(0.5, v / 2);
+ for (int i = 0; i < n; i++) {
+ a = xL + i * dx;
+ b = a + dx;
+ xc = a + dx / 2.0;
+ tmp = 0.0;
+ tmp = tmp + densityfunc(a, v, B) / 3;
+ tmp = tmp + densityfunc(xc, v, B) * 4 / 3;
+ tmp = tmp + densityfunc(b, v, B) / 3;
+ tmp = tmp * dx / 2;
+ goukei += tmp;
+ }
+ return goukei;
+ }
+
+ private static double densityfunc(double x, double v, double B) {
+ return Math.pow(1.0 + (x * x / v), -(v + 1) / 2) / (Math.sqrt(v) * B);
+ }
+
+ /**
+ * うまくいかないので使わない
+ * @param u
+ * @param v
+ * @return
+ */
+ public static double simpson_betafunc(double u, double v) {
+ double xL = 0.0;
+ double xR = 1.0;
+ int n = 100;
+ double goukei = 0.0;
+ double a, b, xc, tmp;
+ double dx = (xR - xL) / n;
+ for (int i = 0; i < n; i++) {
+ a = xL + i * dx;
+ b = a + dx;
+ xc = a + dx / 2.0;
+ tmp = 0.0;
+ tmp = tmp + bfunc(a, u, v) / 3;
+ tmp = tmp + bfunc(xc, u, v) * 4 / 3;
+ tmp = tmp + bfunc(b, u, v) / 3;
+ tmp = tmp * dx / 2;
+ goukei += tmp;
+ System.out.println(i+" "+tmp+" "+goukei);
+ }
+ return goukei;
+ }
+
+ /**
+ * うまくいかないので使わない
+ * @param x
+ * @param u
+ * @param v
+ * @return
+ */
+ private static double bfunc(double x, double u, double v) {
+ double d = Math.pow(x, u - 1) * Math.pow(1.0 - x, v - 1);
+ System.out.println("bfunc "+x+" "+u+" "+v+" = "+d);
+ return d;
+ }
+
}