#include <math.h>
#include "lib/mlrstat.h"

// ================================================================
// These are intended for streaming (i.e. single-pass) applications. Otherwise
// the formulas look different (and are more intuitive).
// ================================================================

// ----------------------------------------------------------------
// Univariate linear regression
// ----------------------------------------------------------------
// There are N (xi, yi) pairs.
//
// minimize E = sum (yi - m xi - b)^2
//
// Set the two partial derivatives to zero and solve for m and b:
//
// DE/Dm = sum 2 (yi - m xi - b) (-xi) = 0
// DE/Db = sum 2 (yi - m xi - b) (-1)  = 0
//
// sum (yi - m xi - b) (xi) = 0
// sum (yi - m xi - b)      = 0
//
// sum (xi yi - m xi^2 - b xi) = 0
// sum (yi - m xi - b)         = 0
//
// m sum(xi^2) + b sum(xi) = sum(xi yi)
// m sum(xi)   + b N       = sum(yi)
//
// [ sum(xi^2)   sum(xi) ] [ m ] = [ sum(xi yi) ]
// [ sum(xi)     N       ] [ b ] = [ sum(yi)    ]
//
// [ m ] = [ sum(xi^2) sum(xi) ]^-1  [ sum(xi yi) ]
// [ b ]   [ sum(xi)   N       ]     [ sum(yi)    ]
//
//       = [ N         -sum(xi)  ]  [ sum(xi yi) ] * 1/D
//         [ -sum(xi)   sum(xi^2)]  [ sum(yi)    ]
//
// where
//
//   D = N sum(xi^2) - sum(xi)^2.
//
// So
//
//      N sum(xi yi) - sum(xi) sum(yi)
// m = --------------------------------
//                   D
//
//      -sum(xi)sum(xi yi) + sum(xi^2) sum(yi)
// b = ----------------------------------------
//                   D
//
// ----------------------------------------------------------------

void mlr_get_linear_regression_ols(unsigned long long n, double sumx, double sumx2, double sumxy, double sumy,
	double* pm, double* pb)
{
	double D =  n * sumx2 - sumx*sumx;
	double m = (n * sumxy - sumx * sumy) / D;
	double b = (-sumx * sumxy + sumx2 * sumy) / D;

	*pm = m;
	*pb = b;
}

// We would need a second pass through the data to compute the error-bars given
// the data and the m and the b.
//
//	# Young 1962, pp. 122-124.  Compute sample variance of linear
//	# approximations, then variances of m and b.
//	var_z = 0.0
//	for i in range(0, N):
//		var_z += (m * xs[i] + b - ys[i])**2
//	var_z /= N
//
//	var_m = (N * var_z) / D
//	var_b = (var_z * sumx2) / D
//
//	output = [m, b, math.sqrt(var_m), math.sqrt(var_b)]

// ----------------------------------------------------------------
double mlr_get_var(unsigned long long n, double sum, double sum2) {
	double mean = sum / n;
	double numerator = sum2 - mean*(2.0*sum - n*mean);
	if (numerator < 0.0) // round-off error
		numerator = 0.0;
	double denominator = n - 1LL;
	return numerator / denominator;
}

// ----------------------------------------------------------------
// Unbiased estimator:
//    (1/n)   sum{(xi-mean)**3}
//  -----------------------------
// [(1/(n-1)) sum{(xi-mean)**2}]**1.5

// mean = sumx / n; n mean = sumx

// sum{(xi-mean)^3}
//   = sum{xi^3 - 3 mean xi^2 + 3 mean^2 xi - mean^3}
//   = sum{xi^3} - 3 mean sum{xi^2} + 3 mean^2 sum{xi} - n mean^3
//   = sumx3 - 3 mean sumx2 + 3 mean^2 sumx - n mean^3
//   = sumx3 - 3 mean sumx2 + 3n mean^3 - n mean^3
//   = sumx3 - 3 mean sumx2 + 2n mean^3
//   = sumx3 - mean*(3 sumx2 + 2n mean^2)

// sum{(xi-mean)^2}
//   = sum{xi^2 - 2 mean xi + mean^2}
//   = sum{xi^2} - 2 mean sum{xi} + n mean^2
//   = sumx2 - 2 mean sumx + n mean^2
//   = sumx2 - 2 n mean^2 + n mean^2
//   = sumx2 - n mean^2

double mlr_get_skewness(unsigned long long n, double sumx, double sumx2, double sumx3) {
	double mean = sumx / n;
	double numerator = sumx3 - mean*(3*sumx2 - 2*n*mean*mean);
	numerator = numerator / n;
	double denominator = (sumx2 - n*mean*mean) / (n-1);
	denominator = pow(denominator, 1.5);
	return numerator / denominator;
}

// Unbiased:
//  (1/n) sum{(x-mean)**4}
//  ----------------------- - 3
// [(1/n) sum{(x-mean)**2}]**2

// sum{(xi-mean)^4}
//   = sum{xi^4 - 4 mean xi^3 + 6 mean^2 xi^2 - 4 mean^3 xi + mean^4}
//   = sum{xi^4} - 4 mean sum{xi^3} + 6 mean^2 sum{xi^2} - 4 mean^3 sum{xi} + n mean^4
//   = sum{xi^4} - 4 mean sum{xi^3} + 6 mean^2 sum{xi^2} - 4 n mean^4 + n mean^4
//   = sum{xi^4} - 4 mean sum{xi^3} + 6 mean^2 sum{xi^2} - 3 n mean^4
//   = sum{xi^4} - mean*(4 sum{xi^3} - 6 mean sum{xi^2} + 3 n mean^3)
//   = sumx4 - mean*(4 sumx3 - 6 mean sumx2 + 3 n mean^3)
//   = sumx4 - mean*(4 sumx3 - mean*(6 sumx2 - 3 n mean^2))

double mlr_get_kurtosis(unsigned long long n, double sumx, double sumx2, double sumx3, double sumx4) {
	double mean = sumx / n;
	double numerator = sumx4 - mean*(4*sumx3 - mean*(6*sumx2 - 3*n*mean*mean));
	numerator = numerator / n;
	double denominator = (sumx2 - n*mean*mean) / n;
	denominator = denominator * denominator;
	return numerator / denominator - 3.0;
}

// ----------------------------------------------------------------
// Non-streaming implementation:
//
// def find_sample_covariance(xs, ys):
//      n = len(xs)
//      mean_x = find_mean(xs)
//      mean_y = find_mean(ys)
//
//      sum = 0.0
//      for k in range(0, n):
//              sum += (xs[k] - mean_x) * (ys[k] - mean_y)
//
//      return sum / (n-1.0)

double mlr_get_cov(unsigned long long n, double sumx, double sumy, double sumxy) {
	double meanx = sumx / n;
	double meany = sumy / n;
	double numerator = sumxy - meanx*sumy - meany*sumx + n*meanx*meany;
	double denominator = n - 1;
	return numerator / denominator;
}

// ----------------------------------------------------------------
void mlr_get_cov_matrix(unsigned long long n,
	double sumx, double sumx2, double sumy, double sumy2, double sumxy,
	double Q[2][2])
{
	double denominator = n - 1;
	Q[0][0] = (sumx2 - sumx*sumx/n) / denominator;
	Q[0][1] = (sumxy - sumx*sumy/n) / denominator;
	Q[1][0] = Q[0][1];
	Q[1][1] = (sumy2 - sumy*sumy/n) / denominator;
}

// ----------------------------------------------------------------
// Principal component analysis can be used for linear regression:
//
// * Compute the covariance matrix for the x's and y's.
//
// * Find its eigenvalues and eigenvectors of the cov. (This is real-symmetric
//   so Jacobi iteration is simple and fine.)
//
// * The principal eigenvector points in the direction of the fit.
//
// * The covariance matrix is computed on zero-mean data so the intercept
//   is zero. The fit equation is of the form (y - nu) = m*(x - mu) where mu
//   and nu are x and y means, respectively.
//
// * If the fit is perfect then the 2nd eigenvalue will be zero; if the fit is
//   good then the 2nd eigenvalue will be smaller; if the fit is bad then
//   they'll be about the same. I use 1 - |lambda2|/|lambda1| as an indication
//   of quality of the fit.
//
// Standard ("ordinary least-squares") linear regression is appropriate when
// the errors are thought to be all in the y's. PCA ("total least-squares") is
// appropriate when the x's and the y's are thought to both have errors.

void mlr_get_linear_regression_pca(
	// Inputs:
	double eigenvalue_1,
	double eigenvalue_2,
	double eigenvector_1[2],
	double eigenvector_2[2],
	double x_mean, double y_mean,
	// Outputs:
	double* pm, double* pb, double* pquality)
{
	double abs_1 = fabs(eigenvalue_1);
	double abs_2 = fabs(eigenvalue_2);
	double quality = 1.0;
	if (abs_1 == 0.0)
		quality = 0.0;
	else if (abs_2 > 0.0)
		quality = 1.0 - abs_2 / abs_1;

	double a0 = eigenvector_1[0];
	double a1 = eigenvector_1[1];
	double m = a1 / a0;
	double b = y_mean - m * x_mean;

	*pm = m;
	*pb = b;
	*pquality = quality;
}