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 = "+p1);
System.out.println("----");
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, 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));
// }
}
/**
* 合計
*/
public static double sum(double[] ary) {
double tmp = 0;
for (int i = 0; i < ary.length; i++) {
tmp += ary[i];
}
return tmp;
}
/**
* 平均
*/
public static double avg(double[] ary) {
double sum = sum(ary);
return sum / ary.length;
}
/**
* 偏差：平均値との差の2乗
*/
public static double hensa(double[] ary) {
double avg = avg(ary);
double tmp = 0;
for (int i = 0; i < ary.length; i++) {
tmp += (ary[i] - avg) * (ary[i] - avg);
}
return tmp;
}
/**
* 不偏分散： 偏差を、データ数n-1 で割ったもの vars
*/
public static double var(double[] ary) {
double hensa = hensa(ary);
return hensa / (ary.length - 1);
}
/**
* 標本分散： 偏差を、データ数n で割ったもの
*/
public static double varp(double[] ary) {
double hensa = hensa(ary);
return hensa / (ary.length);
}
/**
* 差の配列
*
* @param ary1
* @param ary2
* @return
*/
public static double[] diffary(double[] ary1, double[] ary2) {
double[] ret = new double[ary1.length];
for (int i = 0; i < ary1.length; i++) {
ret[i] = ary1[i] - ary2[i];
}
return ret;
}
/**
* 対応のあるt検定の t値
*/
public static double paired_t(double[] ary1, double[] ary2) {
// 差の平均値 / sqrt(差の不偏分散 / データ数)
double avg_sa = avg(diffary(ary1, ary2));
double var_sa = var(diffary(ary1, ary2));
return avg_sa / Math.sqrt(var_sa / ary1.length);
}
/**
* 対応がない(なくてもよい) t検定 Welchのt
*/
public static double welch_t(double[] ary1, double[] ary2) {
// 不偏分散/データ数
double stddev1 = var(ary1) / ary1.length;
double stddev2 = var(ary2) / ary2.length;
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) {
// 不偏分散/データ数
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;
}
/**
* TDIST by commons-math3 TDistribution
*
* @param x
* @param df
* @param tails
* @return
*/
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
*
*
* Calculates the probability from -INF to X under Student's Distribution
* Ported to PHP from Javascript implementation found at
* http://www.math.ucla.edu/~tom/distributions/tDist.html
* // w w w .j av a 2 s. c o m
*
* @param x
* @param deg_free
* @return
*/
public static double studentsCdf(double x, double deg_free) {
if (deg_free <= 0) {
throw new IllegalArgumentException("The degrees of freedom need to be positive.");
}
double A = deg_free / 2.0;
double S = A + 0.5;
double Z = deg_free / (deg_free + x * x);
double BT = Math.exp(logGamma(S) - logGamma(0.5) - logGamma(A) + A * Math.log(Z) + 0.5 * Math.log(1.0 - Z));
double betacdf;
if (Z < (A + 1.0) / (S + 2.0)) {
betacdf = BT * betinc(Z, A, 0.5);
} else {
betacdf = 1 - BT * betinc(1.0 - Z, 0.5, A);
}
double tcdf;
if (x < 0) {
tcdf = betacdf / 2.0;
} else {
tcdf = 1.0 - betacdf / 2.0;
}
return tcdf;
}
/**
* Log Gamma Function
*
* @param Z
* @return
*/
public static double logGamma(double Z) {
double S = 1.0 + 76.18009173 / Z - 86.50532033 / (Z + 1.0) + 24.01409822 / (Z + 2.0)
- 1.231739516 / (Z + 3.0) + 0.00120858003 / (Z + 4.0) - 0.00000536382 / (Z + 5.0);
double LG = (Z - 0.5) * Math.log(Z + 4.5) - (Z + 4.5) + Math.log(S * 2.50662827465);
return LG;
}
/**
* Internal function used by StudentCdf
*
* @param x
* @param A
* @param B
* @return
*/
protected static double betinc(double x, double A, double B) {
double A0 = 0.0;
double B0 = 1.0;
double A1 = 1.0;
double B1 = 1.0;
double M9 = 0.0;
double A2 = 0.0;
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;
B0 = B1 + C9 * B0;
M9 = M9 + 1;
C9 = M9 * (B - M9) * x / (A + 2.0 * M9 - 1.0) / (A + 2.0 * M9);
A1 = A0 + C9 * A1;
B1 = B0 + C9 * B1;
A0 = A0 / B1;
B0 = B0 / B1;
A1 = A1 / B1;
B1 = 1.0;
}
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;
}
}
