xref: /zscilib-3.4.0/src/statistics.c

/*
 * Copyright (c) 2021 Marti Riba Pons
 *
 * SPDX-License-Identifier: Apache-2.0
 */

#include <errno.h>
#include <stdio.h>
#include <stdbool.h>
#include <string.h>
#include <zsl/zsl.h>
#include <zsl/statistics.h>

int zsl_sta_mean(struct zsl_vec *v, zsl_real_t *m)
{
	zsl_vec_ar_mean(v, m);

	return 0;
}

int zsl_sta_trim_mean(struct zsl_vec *v, zsl_real_t p, zsl_real_t *m)
{
#if CONFIG_ZSL_BOUNDS_CHECKS
	/* Make sure p is between 0 and 50. */
	if (p > 50.0 || p < 0.0) {
		return -EINVAL;
	}
#endif

	ZSL_VECTOR_DEF(w, v->sz);
	zsl_real_t per_l, per_h;
	size_t first_val = 0, count = 0;

	*m = 0.0;
	zsl_sta_percentile(v, p, &per_l);
	zsl_sta_percentile(v, 100 - p, &per_h);

	zsl_vec_init(&w);
	zsl_vec_sort(v, &w);

	for (size_t i = 0; i < v->sz; i++) {
		if (w.data[i] >= per_l && w.data[i] <= per_h) {
			count++;
		}
		if (count == 1) {
			first_val = i;
		}
	}

	ZSL_VECTOR_DEF(sub, count);
	zsl_vec_get_subset(&w, first_val, count, &sub);
	zsl_sta_mean(&sub, m);

	return 0;
}

int zsl_sta_weighted_mean(struct zsl_vec *v, struct zsl_vec *w, zsl_real_t *m)
{
	zsl_real_t sumw = 0.0, sumwx = 0.0;

	for (size_t i = 0; i < v->sz; i++) {
		sumw += w->data[i];
	}

#if CONFIG_ZSL_BOUNDS_CHECKS
	/* Make sure the vectors dimensions match. */
	if (v->sz != w->sz) {
		return -EINVAL;
	}

	/* Make sure the weights are positive or zero. */
	for (size_t i = 0; i < v->sz; i++) {
		if (w->data[i] < 0.0) {
			return -EINVAL;
		}
	}

	/* Make sure that not all of the weights are zero. */
	if (ZSL_ABS(sumw) < 1E-6) {
		return -EINVAL;
	}
#endif

	for (size_t i = 0; i < v->sz; i++) {
		sumwx += w->data[i] * v->data[i];
	}

	*m = sumwx / sumw;

	return 0;
}

int zsl_sta_time_weighted_mean(struct zsl_vec *v, struct zsl_vec *t,
			       zsl_real_t *m)
{
#if CONFIG_ZSL_BOUNDS_CHECKS
	/* Make sure the vectors dimensions match. */
	if (v->sz != t->sz) {
		return -EINVAL;
	}
	/* Make sure that the values in 'v' are positive. */
	for (size_t i = 0; i < v->sz; i++) {
		if (v->data[i] < 0.0) {
			return -EINVAL;
		}
	}
	/* The vector 'x' can't have any repeated values. */
	for (size_t i = 0; i < t->sz; i++) {
		if (zsl_vec_contains(t, t->data[i], 1E-6) > 1) {
			return -EINVAL;
		}
	}
#endif

	ZSL_VECTOR_DEF(ts, t->sz);
	ZSL_VECTOR_DEF(vs, v->sz);

	zsl_real_t sum = 0.0;

	zsl_vec_sort(t, &ts);

	for (size_t i = 0; i < t->sz; i++) {
		for (size_t g = 0; g < t->sz; g++) {
			if (ZSL_ABS(ts.data[i] - t->data[g]) < 1E-6) {
				vs.data[i] = v->data[g];
			}
		}
	}

	for (size_t j = 1; j < t->sz; j++) {
		sum += (ts.data[j] - ts.data[j - 1]) *
		       (vs.data[j] + vs.data[j - 1]) / 2.0;
	}

	*m = sum / (ts.data[(t->sz - 1)] - ts.data[0]);

	return 0;
}

int zsl_sta_demean(struct zsl_vec *v, struct zsl_vec *w)
{
#if CONFIG_ZSL_BOUNDS_CHECKS
	/* Make sure the vectors dimensions match. */
	if (v->sz != w->sz) {
		return -EINVAL;
	}
#endif

	zsl_real_t m;

	zsl_vec_copy(w, v);
	zsl_sta_mean(v, &m);
	zsl_vec_scalar_add(w, -m);

	return 0;
}

int zsl_sta_percentile(struct zsl_vec *v, zsl_real_t p, zsl_real_t *val)
{
#if CONFIG_ZSL_BOUNDS_CHECKS
	/* Make sure p is between 0 and 100. */
	if (p > 100.0 || p < 0.0) {
		return -EINVAL;
	}
#endif

	ZSL_VECTOR_DEF(w, v->sz);

	zsl_vec_init(&w);
	zsl_vec_sort(v, &w);

	zsl_real_t x = (p * v->sz) / 100.;
	zsl_real_t per = ZSL_FLOOR(x);

	if (x == per) {
		*val = (w.data[(int)per] + w.data[(int)per - 1]) / 2.;
	} else {
		*val = w.data[(int)per];
	}

	return 0;
}

int zsl_sta_median(struct zsl_vec *v, zsl_real_t *m)
{
	zsl_sta_percentile(v, 50, m);

	return 0;
}

int zsl_sta_weighted_median(struct zsl_vec *v, struct zsl_vec *w, zsl_real_t *m)
{
#if CONFIG_ZSL_BOUNDS_CHECKS
	/* Make sure the vectors dimensions match. */
	if (v->sz != w->sz) {
		return -EINVAL;
	}

	/* Make sure the weights are positive or zero. */
	zsl_real_t sum = 0.0;
	for (size_t i = 0; i < v->sz; i++) {
		sum += w->data[i];
		if (w->data[i] < 0.0) {
			return -EINVAL;
		}
	}

	/* Make sure that that the sum of the weights is 1. */
	if (ZSL_ABS(sum - 1.0) > 1E-6) {
		return -EINVAL;
	}
#endif

	ZSL_VECTOR_DEF(vsort, v->sz);
	zsl_real_t lsum = 0.0;
	size_t i = 0;
	while (lsum <= 0.5) {
		lsum += w->data[i];
		i++;
	}

	lsum -= w->data[i - 1];
	zsl_vec_init(&vsort);
	zsl_vec_sort(v, &vsort);

	if (ZSL_ABS(lsum - 0.5) < 1E-6) {
		*m = (vsort.data[i - 1] + vsort.data[i - 2]) / 2.0;
	} else {
		*m = vsort.data[i - 1];
	}

	return 0;
}

int zsl_sta_quart(struct zsl_vec *v, zsl_real_t *q1, zsl_real_t *q2,
		  zsl_real_t *q3)
{
	zsl_sta_percentile(v, 25, q1);
	zsl_sta_percentile(v, 50, q2);
	zsl_sta_percentile(v, 75, q3);

	return 0;
}

int zsl_sta_quart_range(struct zsl_vec *v, zsl_real_t *r)
{
	zsl_real_t q1, q3;

	zsl_sta_percentile(v, 25, &q1);
	zsl_sta_percentile(v, 75, &q3);

	*r = q3 - q1;

	return 0;
}

int zsl_sta_mode(struct zsl_vec *v, struct zsl_vec *w)
{
#if CONFIG_ZSL_BOUNDS_CHECKS
	/* Make sure w and v are the same length. */
	if (v->sz != w->sz) {
		return -EINVAL;
	}
#endif

	size_t i, j = 0, count, maxcount = 0;

	ZSL_VECTOR_DEF(u, v->sz);
	ZSL_VECTOR_DEF(w_temp, v->sz);

	zsl_vec_init(w);

	for (i = 0; i < v->sz; i++) {
		count = zsl_vec_contains(v, v->data[i], 1E-7);
		u.data[i] = count;
		if (count > maxcount) {
			maxcount = count;
		}
	}

	for (i = 0; i < v->sz; i++) {
		if (u.data[i] == maxcount) {
			w_temp.data[j] = v->data[i];
			j++;
		}
	}

	w_temp.sz = j;
	count = 0;

	for (i = 0; i < j; i++) {
		if (w_temp.data[i] >= 1E-5 || w_temp.data[i] <= 1E-5) {
			if (zsl_vec_contains(w, w_temp.data[i], 1E-5) == 0) {
				w->data[count] = w_temp.data[i];
				count++;
			}
		}
	}

	if (zsl_vec_contains(&w_temp, 0.0, 1E-5) > 0) {
		count++;
	}

	w->sz = count;

	return 0;
}

int zsl_sta_data_range(struct zsl_vec *v, zsl_real_t *r)
{
	ZSL_VECTOR_DEF(w, v->sz);

	zsl_vec_sort(v, &w);

	*r = w.data[v->sz - 1] - w.data[0];

	return 0;
}

int zsl_sta_mean_abs_dev(struct zsl_vec *v, zsl_real_t *m)
{
#if CONFIG_ZSL_BOUNDS_CHECKS
	/* Make sure v has at least dimension 1. */
	if (v->sz == 0) {
		return -EINVAL;
	}
#endif
	ZSL_VECTOR_DEF(vdm, v->sz);
	zsl_sta_demean(v, &vdm);

	zsl_real_t sum = 0.0;
	for (size_t i = 0; i < v->sz; i++) {
		sum += ZSL_ABS(vdm.data[i]);
	}

	*m = sum / v->sz;

	return 0;
}

int zsl_sta_median_abs_dev(struct zsl_vec *v, zsl_real_t *m)
{
	ZSL_VECTOR_DEF(med, v->sz);
	zsl_real_t median;

	zsl_sta_median(v, &median);
	for (size_t i = 0; i < v->sz; i++) {
		med.data[i] = ZSL_ABS(v->data[i] - median);
	}

	zsl_sta_median(&med, m);

	return 0;
}

int zsl_sta_var(struct zsl_vec *v, zsl_real_t *var)
{
	ZSL_VECTOR_DEF(w, v->sz);
	*var = 0;

	zsl_sta_demean(v, &w);

	for (size_t i = 0; i < v->sz; i++) {
		*var += w.data[i] * w.data[i];
	}

	*var /= v->sz - 1;

	return 0;
}

int zsl_sta_std_dev(struct zsl_vec *v, zsl_real_t *s)
{
	zsl_real_t var;

	zsl_sta_var(v, &var);

	*s = ZSL_SQRT(var);

	return 0;
}

int zsl_sta_covar(struct zsl_vec *v, struct zsl_vec *w, zsl_real_t *c)
{
#if CONFIG_ZSL_BOUNDS_CHECKS
	/* Make sure v and w are equal length. */
	if (v->sz != w->sz) {
		return -EINVAL;
	}
#endif

	ZSL_VECTOR_DEF(v_dm, v->sz);
	ZSL_VECTOR_DEF(w_dm, w->sz);

	zsl_sta_demean(v, &v_dm);
	zsl_sta_demean(w, &w_dm);

	*c = 0;
	for (size_t i = 0; i < v->sz; i++) {
		*c += v_dm.data[i] * w_dm.data[i];
	}

	*c /= v->sz - 1;

	return 0;
}

int zsl_sta_covar_mtx(struct zsl_mtx *m, struct zsl_mtx *mc)
{
#if CONFIG_ZSL_BOUNDS_CHECKS
	/* Make sure 'mc' is a square matrix with same num. columns as 'm'. */
	if (mc->sz_rows != mc->sz_cols || mc->sz_cols != m->sz_cols) {
		return -EINVAL;
	}
#endif

	zsl_real_t c;

	ZSL_VECTOR_DEF(v1, m->sz_rows);
	ZSL_VECTOR_DEF(v2, m->sz_rows);

	for (size_t i = 0; i < m->sz_cols; i++) {
		for (size_t j = 0; j < m->sz_cols; j++) {
			zsl_mtx_get_col(m, i, v1.data);
			zsl_mtx_get_col(m, j, v2.data);
			zsl_sta_covar(&v1, &v2, &c);
			mc->data[(i * m->sz_cols) + j] = c;
		}
	}

	return 0;
}

int zsl_sta_linear_reg(struct zsl_vec *x, struct zsl_vec *y,
		       struct zsl_sta_linreg *c)
{
#if CONFIG_ZSL_BOUNDS_CHECKS
	/* Make sure x and y are equal length. */
	if (x->sz != y->sz) {
		return -EINVAL;
	}
#endif

	zsl_real_t sumx, sumy, sumxy, sumxx, sumyy;

	sumx = sumy = sumxy = sumxx = sumyy = 0.0;

	for (size_t i = 0; i < x->sz; i++) {
		sumx += x->data[i];
		sumy += y->data[i];
		sumxy += x->data[i] * y->data[i];
		sumxx += x->data[i] * x->data[i];
		sumyy += y->data[i] * y->data[i];
	}

	c->slope = (x->sz * sumxy - sumx * sumy) / (x->sz * sumxx - sumx * sumx);
	c->intercept = (sumy - c->slope * sumx) / x->sz;
	c->correlation = (x->sz * sumxy - sumx * sumy) /
			 ZSL_SQRT((x->sz * sumxx - sumx * sumx) *
				  (x->sz * sumyy - sumy * sumy));

	return 0;
}

#ifndef CONFIG_ZSL_SINGLE_PRECISION
int zsl_sta_mult_linear_reg(struct zsl_mtx *x, struct zsl_vec *y,
			    struct zsl_vec *b, zsl_real_t *r)
{

#if CONFIG_ZSL_BOUNDS_CHECKS
	/* Make sure matrices and vectors' sizes match. */
	if (x->sz_rows != y->sz || x->sz_cols + 1 != b->sz) {
		return -EINVAL;
	}
#endif

	ZSL_MATRIX_DEF(x_exp, x->sz_rows, (x->sz_cols + 1));
	ZSL_MATRIX_DEF(x_trans, (x->sz_cols + 1), x->sz_rows);
	ZSL_MATRIX_DEF(xx, (x->sz_cols + 1), (x->sz_cols + 1));
	ZSL_MATRIX_DEF(xinv, (x->sz_cols + 1), (x->sz_cols + 1));
	ZSL_MATRIX_DEF(ymtx, y->sz, 1);
	ZSL_MATRIX_DEF(xtemp, (x->sz_cols + 1), x->sz_rows);
	ZSL_MATRIX_DEF(bmtx, (x->sz_cols + 1), 1);
	ZSL_MATRIX_DEF(xtemp2, x->sz_rows, 1);
	ZSL_MATRIX_DEF(emtx, x->sz_rows, 1);

	zsl_real_t v[x->sz_rows];
	zsl_mtx_init(&x_exp, NULL);
	for (size_t i = 0; i < x->sz_rows; i++) {
		v[i] = 1.;
	}

	zsl_mtx_set_col(&x_exp, 0, v);
	for (size_t j = 0; j < x->sz_cols; j++) {
		zsl_mtx_get_col(x, j, v);
		zsl_mtx_set_col(&x_exp, j + 1, v);
	}

	zsl_mtx_trans(&x_exp, &x_trans);
	zsl_mtx_mult(&x_trans, &x_exp, &xx);

	zsl_real_t det;
	zsl_mtx_deter(&xx, &det);
	if (ZSL_ABS(det) < 1E-6) {
		/*
		 * Currently limited to square matrices. pinv could be used,
		 * but is too resource-intensive to add at the momemt.
		 */
		return -EINVAL;
	} else {
		zsl_mtx_inv(&xx, &xinv);
	}

	zsl_mtx_from_arr(&ymtx, y->data);
	zsl_mtx_mult(&xinv, &x_trans, &xtemp);
	zsl_mtx_mult(&xtemp, &ymtx, &bmtx);
	zsl_vec_from_arr(b, bmtx.data);

	zsl_mtx_mult(&x_exp, &bmtx, &xtemp2);
	zsl_mtx_sub(&ymtx, &xtemp2, &emtx);
	zsl_real_t e_norm = 0.0, ymean, ysum = 0.0;
	zsl_sta_mean(y, &ymean);
	for (size_t g = 0; g < x->sz_rows; g++) {
		e_norm += emtx.data[g] * emtx.data[g];
		ysum += (y->data[g] - ymean) * (y->data[g] - ymean);
	}

	*r = 1. - e_norm / ysum;

	return 0;
}
#endif

#ifndef CONFIG_ZSL_SINGLE_PRECISION
int zsl_sta_weighted_mult_linear_reg(struct zsl_mtx *x, struct zsl_vec *y,
				     struct zsl_vec *w, struct zsl_vec *b, zsl_real_t *r)
{
#if CONFIG_ZSL_BOUNDS_CHECKS
	/* Make sure matrices and vectors' sizes match. */
	if (x->sz_rows != y->sz || x->sz_cols + 1 != b->sz || x->sz_rows != w->sz) {
		return -EINVAL;
	}
	/* Make sure no weight is zero. */
	for (size_t k = 0; k < w->sz; k++) {
		if (ZSL_ABS(w->data[k]) < 1E-6) {
			return -EINVAL;
		}
	}
#endif

	ZSL_MATRIX_DEF(x_exp, x->sz_rows, (x->sz_cols + 1));
	ZSL_MATRIX_DEF(x_trans, (x->sz_cols + 1), x->sz_rows);
	ZSL_MATRIX_DEF(xw, (x->sz_cols + 1), x->sz_rows);
	ZSL_MATRIX_DEF(xx, (x->sz_cols + 1), (x->sz_cols + 1));
	ZSL_MATRIX_DEF(xinv, (x->sz_cols + 1), (x->sz_cols + 1));
	ZSL_MATRIX_DEF(ymtx, y->sz, 1);
	ZSL_MATRIX_DEF(xtemp, (x->sz_cols + 1), x->sz_rows);
	ZSL_MATRIX_DEF(bmtx, (x->sz_cols + 1), 1);
	ZSL_MATRIX_DEF(xtemp2, x->sz_rows, 1);
	ZSL_MATRIX_DEF(emtx, x->sz_rows, 1);
	ZSL_MATRIX_DEF(idx, w->sz, w->sz);

	zsl_mtx_init(&idx, zsl_mtx_entry_fn_identity);
	for (size_t k = 0; k < w->sz; k++) {
		w->data[k] = 1. / w->data[k];
		zsl_mtx_scalar_mult_row_d(&idx, k, w->data[k]);
	}

	zsl_real_t v[x->sz_rows];
	zsl_mtx_init(&x_exp, NULL);
	for (size_t i = 0; i < x->sz_rows; i++) {
		v[i] = 1.;
	}

	zsl_mtx_set_col(&x_exp, 0, v);
	for (size_t j = 0; j < x->sz_cols; j++) {
		zsl_mtx_get_col(x, j, v);
		zsl_mtx_set_col(&x_exp, j + 1, v);
	}

	zsl_mtx_trans(&x_exp, &x_trans);
	zsl_mtx_mult(&x_trans, &idx, &xw);
	zsl_mtx_mult(&xw, &x_exp, &xx);

	zsl_real_t det;
	zsl_mtx_deter(&xx, &det);
	if (ZSL_ABS(det) < 1E-6) {
		/*
		 * Currently limited to square matrices. pinv could be used,
		 * but is too resource-intensive to add at the momemt.
		 */
		return -EINVAL;
	} else {
		zsl_mtx_inv(&xx, &xinv);
	}

	zsl_mtx_from_arr(&ymtx, y->data);
	zsl_mtx_mult(&xinv, &xw, &xtemp);
	zsl_mtx_mult(&xtemp, &ymtx, &bmtx);
	zsl_vec_from_arr(b, bmtx.data);

	zsl_mtx_mult(&x_exp, &bmtx, &xtemp2);
	zsl_mtx_sub(&ymtx, &xtemp2, &emtx);
	zsl_real_t e_norm = 0.0, ymean, ysum = 0.0;
	zsl_sta_mean(y, &ymean);
	for (size_t g = 0; g < x->sz_rows; g++) {
		e_norm += emtx.data[g] * emtx.data[g];
		ysum += (y->data[g] - ymean) * (y->data[g] - ymean);
	}

	*r = 1. - e_norm / ysum;

	return 0;
}
#endif

#ifndef CONFIG_ZSL_SINGLE_PRECISION
int zsl_sta_quad_fit(struct zsl_mtx *m, struct zsl_vec *b)
{
#if CONFIG_ZSL_BOUNDS_CHECKS
	/* Make sure matrices and vectors' sizes are adequate. */
	if (m->sz_cols != 3 || b->sz != 9) {
		return -EINVAL;
	}
#endif

	ZSL_MATRIX_DEF(x, m->sz_rows, 9);
	ZSL_VECTOR_DEF(xv, 9);
	ZSL_VECTOR_DEF(mv, 3);
	ZSL_MATRIX_DEF(y, m->sz_rows, 1);

	for (size_t i = 0; i < m->sz_rows; i++) {
		zsl_mtx_get_row(m, i, mv.data);
		xv.data[0] = mv.data[0] * mv.data[0];
		xv.data[1] = mv.data[1] * mv.data[1];
		xv.data[2] = mv.data[2] * mv.data[2];
		xv.data[3] = 2.0 * mv.data[0] * mv.data[1];
		xv.data[4] = 2.0 * mv.data[0] * mv.data[2];
		xv.data[5] = 2.0 * mv.data[1] * mv.data[2];
		xv.data[6] = 2.0 * mv.data[0];
		xv.data[7] = 2.0 * mv.data[1];
		xv.data[8] = 2.0 * mv.data[2];
		zsl_mtx_set_row(&x, i, xv.data);
		y.data[i] = 1.0;
	}

	ZSL_MATRIX_DEF(xt, 9, m->sz_rows);
	ZSL_MATRIX_DEF(xtx, 9, 9);
	ZSL_MATRIX_DEF(inv, 9, 9);
	ZSL_MATRIX_DEF(xtmp, 9, (m->sz_rows));
	ZSL_MATRIX_DEF(btmp, 9, 1);

	zsl_mtx_trans(&x, &xt);
	zsl_mtx_mult(&xt, &x, &xtx);
	zsl_mtx_inv(&xtx, &inv);
	zsl_mtx_mult(&inv, &xt, &xtmp);
	zsl_mtx_mult(&xtmp, &y, &btmp);

	zsl_vec_from_arr(b, btmp.data);

	return 0;
}
#endif

int zsl_sta_abs_err(zsl_real_t *val, zsl_real_t *exp_val, zsl_real_t *err)
{
	*err = ZSL_ABS(*val - *exp_val);

	return 0;
}

int zsl_sta_rel_err(zsl_real_t *val, zsl_real_t *exp_val, zsl_real_t *err)
{
	*err = ZSL_ABS(100. * *val - 100. * *exp_val) / *exp_val;

	return 0;
}

int zsl_sta_sta_err(struct zsl_vec *v, zsl_real_t *err)
{
#if CONFIG_ZSL_BOUNDS_CHECKS
	/* Make sure v has at least dimension 1. */
	if (v->sz == 0) {
		return -EINVAL;
	}
#endif

	zsl_real_t var;
	zsl_sta_var(v, &var);
	*err = ZSL_SQRT(var / v->sz);

	return 0;
}