var pi = 3.14159265358979324;

// http://www2.toyo.ac.jp/~aizawa-sho/bunsan/bunsan.html
function statistics(g)
{
	var dn = g.length + 0.0;
	
	var total = 0.0;
	for(var i=0; i<dn; i++) {
		total += g[i];
	}
	var ave = total/dn;
	
	var dx2 = 0;
	for(var i=0; i<dn; i++) {
		var dx = g[i] - ave;
		dx2 += (dx * dx);
	}
	var s2 = dx2/(dn-1.0);
	
	return {mean: ave, sigma2: s2, n: dn};
}

// http://kusuri-jouhou.com/statistics/tkentei2.html
function ttest(g1, g2)
{
	var s1 = statistics(g1);
	var s2 = statistics(g2);
	var c = _ttest(s1, s2);
	var r = {};
	r.detail = "p="+sprintf("%12.9f",c.p);		// with sprintf
	//r.detail = "p="+c.p;						// without sprintf
	if (c.p < 0.05) {
		r.result = "different";
	} else {
		r.result = "same";
	}
	return r;
}

function _ttest(s1, s2)
{
	// calcurate t
	var rt = (s1.mean - s2.mean)/(Math.sqrt( (s1.sigma2/s1.n) + (s2.sigma2/s2.n) ));
	
	// calcurate df
	var df1a = (s1.sigma2/s1.n + s2.sigma2/s2.n);
	var df1 = df1a*df1a;
	var df2a = (s1.sigma2/s1.n);
	var df2b = (s2.sigma2/s2.n);
	var df2 = ( (df2a*df2a)/(s1.n-1) + (df2b*df2b)/(s2.n-1) );
	var rdf = df1/df2;
	
	return {t: rt, df: rdf, p: txp(rt, rdf)};
}

// http://aoki2.si.gunma-u.ac.jp/Hanasi/Algo/getprob.html
function gamma(x)
{
	// Lanczos
	var c0 = 1.000000000178;
	var c1 = 76.180091729406;
	var c2 = -86.505320327112;
	var c3 = 24.01409822223;
	var c4 = -1.23173951614;
	var c5 = 0.00120858003;
	var c6 = -0.00000536382;

	var Ax = c0 + c1/x + c2/(x + 1.0) + c3/(x + 2.0) + c4/(x + 3.0) + c5/(x + 4.0) + c6/(x + 5.0);
	return Math.sqrt(2.0 * pi) * Math.pow( x + 4.5, x - 0.5 ) * Math.exp(-x - 4.5) * Ax;
}

function beta(x, y)
{
	return (gamma(x) * gamma(y))/gamma(x + y);
}

function betai(x, a, b)
{
	var bt = 0.0;

	if (x < 0.0 || x > 1.0) {
		console.log("Bad x in routine betai\n");
		return 99999;
	}
	
	if (x == 0.0 || x == 1.0) {
		bt = 0.0;
	} else {
		bt = Math.exp(gammln(a+b)-gammln(a)-gammln(b)+a*Math.log(x)+b*Math.log(1.0-x));
	}
	
	if (x < (a+1.0)/(a+b+2.0)) {
		return bt*betacf(x, a, b)/a;
	}
	else {
		return 1.0-bt*betacf(1.0-x, b, a)/b;
	}
}

function betacf(x, a, b)
{
	var MAXIT = 100;
	var EPS = 3e-7;
	var FPMIN = 1e-30;
	var m, m2;
	var aa, c, d, del, h, qab, qam, qap;

	qab = a+b;
	qap = a+1;
	qam = a-1;
	c = 1;
	d = 1-qab*x/qap;
	if (Math.abs(d) < FPMIN) d = FPMIN;
	d = 1/d;
	h = d;
	for (m = 1; m <= MAXIT; m++) {
		m2 = 2*m;
		aa = m*(b-m)*x/((qam+m2)*(a+m2));
		d = 1+aa*d;
		if (Math.abs(d) < FPMIN) d = FPMIN;
		c = 1+aa/c;
		if (Math.abs(c) < FPMIN) c = FPMIN;
		d = 1/d;
		h *= d*c;
		aa = -(a+m)*(qab+m)*x/((a+m2)*(qap+m2));
		d = 1+aa*d;
		if (Math.abs(d) < FPMIN) d = FPMIN;
		c = 1+aa/c;
		if (Math.abs(c) < FPMIN) c = FPMIN;
		d = 1/d;
		del = d*c;
		h *= del;
		if (Math.abs(del-1) < EPS) break;
	}
	if (m > MAXIT) {
		console.log("a or b too big, or MAXIT too small in betacf\n");
		return 99999;
	} else {
		return h;
	}
}

function gammln(xx)
{
	var x, y, tmp, ser;
	var j;
	var cof = [
		76.18009172947146,
		-86.50532032941677,
		24.01409824083091,
		-1.231739572450155,
		0.1208650973866179e-2,
		-0.5395239384953e-5
	];

	y = x = xx;
	tmp = x+5.5;
	tmp -= (x+0.5)*Math.log(tmp);
	ser = 1.000000000190015;
	for (j = 0; j <= 5; j++) {
		++y;
		ser = ser + cof[j]/y;
	}
	return -tmp+Math.log(2.5066282746310005*ser/x);
}

function gammp(x, a)
{
	if (x < 0.0 || a <= 0.0) {
		console.log("Invalid arguments in routine gammp");
		return 99999;
	} else if (x < (a+1.0)) {
		return gser(a, x);
	} else {
		return 1-gcf(a, x);
	}
}

function gcf(a, x)
{
	var i, an, b, c, d, del, h, gln;
	var MAXIT, EPS, FPMIN;

	MAXIT = 100;
	EPS = 3.0e-7;
	FPMIN = 1.0e-30;

	gln = gammln(a);
	b = x+1.0-a;
	c = 1.0/FPMIN;
	d = 1.0/b;
	h = d;
	for (i = 1; i <= MAXIT; i++) {
		an = -i*(i-a);
		b += 2.0;
		d = an*d+b;
		if (Math.abs(d) < FPMIN) d = FPMIN;
		c = b+an/c;
		if (Math.abs(c) < FPMIN) c = FPMIN;
		d = 1.0/d;
		del = d*c;
		h *= del;
		if (Math.abs(del-1.0) < EPS) break;
	}
	if (i > MAXIT) console.log("a too large, MAXIT too small in gcf");
	return Math.exp(-x+a*Math.log(x)-gln)*h;
}

function gser(a, x)
{
	var n;
	var sum, del, ap, gln;
	var MAXIT, EPS;

	MAXIT = 100;
	EPS = 3.0e-7;

	gln = gammln(a);
	if (x <= 0.0) {
		if (x < 0.0) console.log("x less than 0 in routine gser");
		return 99999;
	} else {
		ap = a;
		del = sum = 1.0/a;
		for (n = 1; n <= MAXIT; n++) {
			++ap;
			del *= x/ap;
			sum += del;
			if (Math.abs(del) < Math.abs(sum)*EPS) {
				return sum*Math.exp(-x+a*Math.log(x)-gln);
			}
		}
		console.log("a too large, MAXIT too small in routine gser");
		return 99999;
	}
}

function txp(t, df)
{
	var i = 0;
	var p = 0;
	var w = 0;
	var m2pi = 0.636619772367581343076;
	var t1 = Math.abs(t)/Math.sqrt(df)
	var t2 = 1.0/(1.0+t1*t1)
	if ((df % 2) == 0) {
		w = t1 * Math.sqrt(t2)
		p = 1.0-w
		for (i = 2; i <= df-2; i += 2) {
			p -= w *= t2*(i-1)/i
		}
	} else {
		p = 1.0 - m2pi * Math.atan2(t1, 1)
		if (df >= 3) {
			w = m2pi * t1 * t2
			p -= w
			for (i = 3; i <= df-2; i += 2) {
				p -= w *= t2*(i-1)/i
			}
		}
	}
	return (p < 0.0 && Math.abs(p) < 1e-10) ? 0.0 : p
}

// upper probability of normal(gaussian) distribution
function gxp(x)
{
	var pi2 = 0.398942280401432677940;
	var is = -1;
	var y = Math.abs(x);
	var c = y*y;
	var p = 0;
	var z = Math.exp(-c*0.5)*pi2;
	var i = 0;
	if (y < 2.5) {
		for(i=20; i>0; i--) {
			p = (i*c)/(i*2+1+is*p);
			is = -is;
		}
		p = 0.5 - (z * y)/(1.0 - p);
	} else {
		for(i=20; i>0; i--) {
			p = i/(y+p);
		}
		p = z/(y+p);
	}
	return (x < 0.0) ? 1.0-p : p;
}

function pxg(p)
{
	var a =  0.161614603871299e-15
	var b = -0.178837229464158e-13
	var c =  0.888943322341436e-12
	var d = -0.261195052190738e-10
	var e =  0.500533655324426e-09
	var f = -0.646627744010232e-08
	var g =  0.553751002988654e-07
	var h = -0.283629174340630e-06
	var i =  0.539877747512783e-06
	var j =  0.122127251950965e-05
	var k =  0.154574365345128e-04
	var l = -0.234768842092152e-03
	var m = -0.830686207721635e-03
	var n =  0.370684581838722e-01
	var o =  1.57079636931376
	var y = -Math.log(4*p*(1-p))
	return Math.sqrt(((((((((((((((y*a+b)*y+c)*y+d)*y+e)*y+f)*y+g)*y+h)*y+i)*y+j)*y+k)*y+l)*y+m)*y+n)*y+o)*y)
}

function fxp(f, df1, df2)
{
	return betai(df2/(df2+df1*f), df2/2, df1/2);
}

// http://kusuri-jouhou.com/statistics/mann.html
// http://kusuri-jouhou.com/statistics/bunpuhyou2.html
function utest_distribution(n1, n2) {
	if (n1 > n2) {
		var t = n1;
		n1 = n2;
		n2 = t;
	}
	
	//console.log("utest_distribution(n1="+n1+",n2="+n2+")");
	d = -1;
	switch(n1) {
	case  5:
		switch(n2) {
		case  5: d =  2; break;
		case 10: d =  8; break;
		}
		break;
	case  6:
		switch(n2) {
		case  6: d =  5; break;
		case  8: d =  8; break;
		}
		break;
	case 10:
		switch(n2) {
		case 10: d = 23; break;
		}
		break;
	}
	return d;
}

function utest(g1, g2)
{
	var u = 0;
	for(var i=0; i<g1.length; i++) {
		var u1 = 0;
		var u2 = 0;
		for(var j=0; j<g2.length; j++) {
			if (g1[i] < g2[j]) {
				u1++;
			}
			if (g1[i] <= g2[j]) {
				u2++;
			}
		}
		u += (u1 + u2)/2.0;
	}
	
	var r = {};
	if ((g1.length <= 20) && (g2.length <= 20)) {
		//console.log("utest1");
		var c = utest_distribution(g1.length, g2.length);
		//console.log("u="+u+",c="+c);
		r.detail = "u="+u+", c="+c;
		if (u > c) {
			r.result = "same";
		} else {
			r.result = "different";
		}
	} else {
		//console.log("utest2");
		var m = (g1.length * g2.length)/2;
		var s = Math.sqrt( ( (g1.length*g2.length)*(g1.length+g2.length+1) ) / 12.0 );
		var z = (u - m) / s;
		p = gxp(Math.abs(z))*2;
		//console.log("u="+u+",m="+m+",s="+s+",z="+z+", p="+p);
		r.detail = "u="+u+", m="+m+",s="+s+",z="+z+", p="+String(p).substring(0,8);
		if (p < 0.05) {
			r.result = "different";
		} else {
			r.result = "same";
		}
	}
	return r;
}

function ftest(x0, x1)
{
	var s0 = statistics(x0);
	var s1 = statistics(x1);
	var F = s1.sigma2/s0.sigma2;
	var p = fxp(F, x0.length-1, x1.length-1);
	var r = {};
	r.detail = "F="+String(F)+",p="+String(p);
	if (p < 0.05) {
		r.result = "different";
	} else {
		r.result = "same";
	}
	return r;
}

function linearRegression(dy)
{
	var n = dy.length;
	var dx = new Array();
	for(var i=0; i<n; i++) {
		dx.push(i);
	}
	
	var sxy = 0;
	var sx2 = 0;
	for(var i=0; i<n; i++) {
		sxy += dx[i]*dy[i];
		sx2 += dx[i]*dx[i];
	}
	
	var sx = statistics(dx);
	var sy = statistics(dy);
	var c1 = (sxy - n * sx.mean * sy.mean) + 0.0;
	var c0 = (sx2 - n * sx.mean * sx.mean) + 0.0;
	var a = c1/c0;
	var b = sy.mean - a * sx.mean;
	return {slope: a, bias: b};
}

// testing stuff

function test_basic()
{
	// test
	// http://www2.toyo.ac.jp/~aizawa-sho/bunsan/bunsan.html
	var tg = [151, 154, 160, 160, 163, 156, 158, 156, 154, 160];
	var r = statistics(tg);
	console.log("mean="+r.mean+", sigma2="+r.sigma2+", n="+r.n);
}

function test_ttest()
{
	// http://hs-www.hyogo-dai.ac.jp/~kawano/HStat/?2010%2F14th%2FWelch%27s_Test
	// different
	var r1 = { mean: 3.67, sigma2: (1.79*1.79), n: 24 };
	var r2 = { mean: 2.77, sigma2: (1.29*1.29), n: 48 };
	r = _ttest(r1, r2);
	console.log("ttest t="+r.t+", df="+r.df+", p="+r.p);
	
	// test
	var g = new Array();
	g[0] = [9.034078, 8.934082, 8.771254, 8.808801, 8.928719];
	g[1] = [9.273499, 9.361234, 9.181932, 9.207985, 9.180017];
	r = ttest(g[0], g[1]);
	console.log("actual t.result="+r.result+", detail="+r.detail);
	
	console.log("gxp(2.14)="+gxp(2.14)*2);
	console.log("pxg(2.14)="+pxg(0.05));
	
	// http://office.microsoft.com/ja-jp/excel-help/HP010335701.aspx
	g[0] = [  3,  4,  5,  8,  9,  1,  2,  4,  5];
	g[1] = [  6, 19,  3,  2, 14,  4,  5, 17,  1];
	r = ttest(g[0], g[1]);
	console.log("actual t.result="+r.result+", detail="+r.detail);
}

function test_utest() {
	var g = new Array();
	
	// http://kusuri-jouhou.com/statistics/mann.html
	// same
	g[0] = [20, 18, 15, 13, 10, 6];
	g[1] = [17, 16, 12,  9,  8, 4, 4, 2];
	r = utest(g[0], g[1]);
	console.log("case1: result="+r.result+", detail="+r.detail);
	
	// different
	g[0] = [50, 49, 49, 47, 45, 45, 43, 43, 43, 41, 38, 37, 36, 35, 34, 33, 33, 31, 31, 30];
	g[1] = [48, 48, 46, 44, 42, 42, 42, 41, 41, 36, 34, 32, 30, 30, 29, 28, 28, 28, 26, 25, 24, 22];
	r = utest(g[0], g[1]);
	console.log("case2: result="+r.result+", detail="+r.detail);
	
	// test
	g[0] = [9.034078, 8.934082, 8.771254, 8.808801, 8.928719];
	g[1] = [9.273499, 9.361234, 9.181932, 9.207985, 9.180017];
	r = utest(g[0], g[1]);
	console.log("actual: result="+r.result+", detail="+r.detail);
}

function test_ftest() {
	var g = new Array();
	
	g[0] = [49, 40, 52, 37, 55, 38, 45];
	g[1] = [60, 52, 68, 55, 65, 47, 45, 62, 53];
	r = ftest(g[0], g[1]);
	console.log("r.detail="+r.detail);
	console.log("r.result="+r.result);
}

//test_ftest();
