stat.c

Go to the documentation of this file.

/*
 * This work is licensed under a Creative Commons 
 * Attribution-Noncommercial-Share Alike 3.0 United States License.
 * 
 *    http://creativecommons.org/licenses/by-nc-sa/3.0/us/
 * 
 * You are free:
 * 
 *    to Share - to copy, distribute, display, and perform the work
 *    to Remix - to make derivative works
 * 
 * Under the following conditions:
 * 
 *    Attribution. You must attribute the work in the manner specified by the
 *    author or licensor (but not in any way that suggests that they endorse you
 *    or your use of the work).
 * 
 *    Noncommercial. You may not use this work for commercial purposes.
 * 
 *    Share Alike. If you alter, transform, or build upon this work, you may
 *    distribute the resulting work only under the same or similar license to
 *    this one.
 * 
 * For any reuse or distribution, you must make clear to others the license
 * terms of this work. The best way to do this is by including this header.
 * 
 * Any of the above conditions can be waived if you get permission from the
 * copyright holder.
 * 
 * Apart from the remix rights granted under this license, nothing in this
 * license impairs or restricts the author's moral rights.
 */


#include <jwsc/config.h>

#include <stdlib.h>
#include <math.h>
#include <inttypes.h>

#include "jwsc/base/error.h"
#include "jwsc/vector/vector.h"
#include "jwsc/vector/vector_math.h"
#include "jwsc/matrix/matrix.h"
#include "jwsc/matrix/matrix_math.h"
#include "jwsc/stat/stat.h"


void sample_mean_f(float* mean_out, const Vector_f* x)
{
    uint32_t n, N;
    float sum;

    N = x->num_elts;
    sum = 0;

    for (n = 0; n < N; n++)
    {
        sum += x->elts[ n ];
    }

    *mean_out = sum / (float)N;
}

void sample_mean_d(double* mean_out, const Vector_d* x)
{
    uint32_t n, N;
    double sum;

    N = x->num_elts;
    sum = 0;

    for (n = 0; n < N; n++)
    {
        sum += x->elts[ n ];
    }

    *mean_out = sum / (double)N;
}

void sample_variance_f(float* var_out, const Vector_f* x)
{
    uint32_t n, N;
    float sum, ssum;

    N = x->num_elts;
    sum = 0;
    ssum = 0;

    for (n = 0; n < N; n++)
    {
        sum  += x->elts[ n ];
        ssum += x->elts[ n ] * x->elts[ n ];
    }

    if (N > 1)
    {
        *var_out = (1.0 / (float)(N-1)) * (ssum - ((sum*sum) / (float)N));
    }
    else
    {
        *var_out = ssum - sum*sum;
    }
}

void sample_variance_d(double* var_out, const Vector_d* x)
{
    uint32_t n, N;
    double sum, ssum;

    N = x->num_elts;
    sum = 0;
    ssum = 0;

    for (n = 0; n < N; n++)
    {
        sum  += x->elts[ n ];
        ssum += x->elts[ n ] * x->elts[ n ];
    }

    if (N > 1)
    {
        *var_out = (1.0 / (double)(N-1)) * (ssum - ((sum*sum) / (double)N));
    }
    else
    {
        *var_out = ssum - sum*sum;
    }
}

void sample_error_f(float* err_out, const Vector_f* x)
{
    float var;

    sample_variance_f(&var, x);

    *err_out = sqrt(var / (float)(x->num_elts));
}

void sample_error_d(double* err_out, const Vector_d* x)
{
    double var;

    sample_variance_d(&var, x);

    *err_out = sqrt(var / (double)(x->num_elts));
}

void sample_stats_f
(
    float*          mean_out, 
    float*          var_out,
    float*          err_out, 
    const Vector_f* x
)
{
    uint32_t n, N;
    float sum, ssum;

    N = x->num_elts;
    sum = 0;
    ssum = 0;

    for (n = 0; n < N; n++)
    {
        sum  += x->elts[ n ];
        ssum += x->elts[ n ] * x->elts[ n ];
    }

    *mean_out = sum / (float)N;

    if (N > 1)
    {
        *var_out = (1.0 / (float)(N-1)) * (ssum - ((sum*sum) / (float)N));
    }
    else
    {
        *var_out = ssum - sum*sum;
    }

    *err_out = sqrt(*var_out / (float)N);
}

void sample_stats_d
(
    double*         mean_out, 
    double*         var_out, 
    double*         err_out, 
    const Vector_d* x
)
{
    uint32_t n, N;
    double sum, ssum;

    N = x->num_elts;
    sum = 0;
    ssum = 0;

    for (n = 0; n < N; n++)
    {
        sum  += x->elts[ n ];
        ssum += x->elts[ n ] * x->elts[ n ];
    }

    *mean_out = sum / (double)N;

    if (N > 1)
    {
        *var_out = (1.0 / (double)(N-1)) * (ssum - ((sum*sum) / (double)N));
    }
    else
    {
        *var_out = ssum - sum*sum;
    }

    *err_out = sqrt(*var_out / (double)N);
}

void ind_mv_sample_stats_f
(
    Vector_f**      means_out, 
    Vector_f**      vars_out,
    Vector_f**      errs_out, 
    const Matrix_f* x
)
{
    uint32_t d, D;
    uint32_t n, N;
    float sum, ssum;

    D = x->num_cols;
    N = x->num_rows;

    create_vector_f(means_out, D);
    create_vector_f(vars_out, D);
    create_vector_f(errs_out, D);

    for (d = 0; d < D; d++)
    {
        sum = 0;
        ssum = 0;

        for (n = 0; n < N; n++)
        {
            sum  += x->elts[ n ][ d ];
            ssum += x->elts[ n ][ d ] * x->elts[ n ][ d ];
        }

        (*means_out)->elts[ d ] = sum / (float)N;

        if (N > 1)
        {
            (*vars_out)->elts[ d ] = (1.0 / (float)(N-1)) * 
                (ssum - ((sum*sum) / (float)N));
        }
        else
        {
            (*vars_out)->elts[ d ] = ssum - sum*sum;
        }

        (*errs_out)->elts[ d ] = sqrt((*vars_out)->elts[ d ] / (float)N);
    }
}

void ind_mv_sample_stats_d
(
    Vector_d**      means_out, 
    Vector_d**      vars_out, 
    Vector_d**      errs_out, 
    const Matrix_d* x
)
{
    uint32_t d, D;
    uint32_t n, N;
    double sum, ssum;

    D = x->num_cols;
    N = x->num_rows;

    create_vector_d(means_out, D);
    create_vector_d(vars_out, D);
    create_vector_d(errs_out, D);

    for (d = 0; d < D; d++)
    {
        sum = 0;
        ssum = 0;

        for (n = 0; n < N; n++)
        {
            sum  += x->elts[ n ][ d ];
            ssum += x->elts[ n ][ d ] * x->elts[ n ][ d ];
        }

        (*means_out)->elts[ d ] = sum / (double)N;

        if (N > 1)
        {
            (*vars_out)->elts[ d ] = (1.0 / (double)(N-1)) * 
                (ssum - ((sum*sum) / (double)N));
        }
        else
        {
            (*vars_out)->elts[ d ] = ssum - sum*sum;
        }

        (*errs_out)->elts[ d ] = sqrt((*vars_out)->elts[ d ] / (double)N);
    }
}

void mv_sample_mean_f(Vector_f** mean_out, const Matrix_f* x)
{
    uint32_t n, N;
    uint32_t d, D;
    Vector_f* mean;

    D = x->num_cols;
    N = x->num_rows;

    create_zero_vector_f(mean_out, D);
    mean = *mean_out;

    for (n = 0; n < N; n++)
    {
        for (d = 0; d < D; d++)
        {
            mean->elts[ d ] += x->elts[ n ][ d ];
        }
    }

    multiply_vector_by_scalar_f(&mean, mean, 1.0 / (float)N);
}

void mv_sample_mean_d(Vector_d** mean_out, const Matrix_d* x)
{
    uint32_t n, N;
    uint32_t d, D;
    Vector_d* mean;

    D = x->num_cols;
    N = x->num_rows;

    create_zero_vector_d(mean_out, D);
    mean = *mean_out;

    for (n = 0; n < N; n++)
    {
        for (d = 0; d < D; d++)
        {
            mean->elts[ d ] += x->elts[ n ][ d ];
        }
    }

    multiply_vector_by_scalar_d(&mean, mean, 1.0 / (double)N);
}

void mv_sample_covariance_f(Matrix_f** cov_out, const Matrix_f* x)
{
    uint32_t  n, N;
    uint32_t  d, D;
    Matrix_f* cov  = 0;
    Vector_f* mean = 0;
    Matrix_f* m    = 0;
    Matrix_f* m_t  = 0;
    Matrix_f* M    = 0;

    mv_sample_mean_f(&mean, x);

    D = x->num_cols;
    N = x->num_rows;

    create_zero_matrix_f(cov_out, D, D);
    cov = *cov_out;

    create_matrix_f(&m, D, 1);
    create_matrix_f(&m_t, 1, D);
    create_matrix_f(&M, D, D);

    for (n = 0; n < N; n++)
    {
        for (d = 0; d < D; d++)
        {
            m->elts[d][0] = m_t->elts[0][d] = x->elts[n][d] - mean->elts[d];
        }

        multiply_matrices_f(&M, m, m_t);
        add_matrices_f(&cov, cov, M);
    }

    multiply_matrix_by_scalar_f(&cov, cov, 1.0/(double)(N - 1));

    free_vector_f(mean);
    free_matrix_f(m);
    free_matrix_f(m_t);
    free_matrix_f(M);
}

void mv_sample_covariance_d(Matrix_d** cov_out, const Matrix_d* x)
{
    uint32_t  n, N;
    uint32_t  d, D;
    Matrix_d* cov  = 0;
    Vector_d* mean = 0;
    Matrix_d* m    = 0;
    Matrix_d* m_t  = 0;
    Matrix_d* M    = 0;

    mv_sample_mean_d(&mean, x);

    D = x->num_cols;
    N = x->num_rows;

    create_zero_matrix_d(cov_out, D, D);
    cov = *cov_out;

    create_matrix_d(&m, D, 1);
    create_matrix_d(&m_t, 1, D);
    create_matrix_d(&M, D, D);

    for (n = 0; n < N; n++)
    {
        for (d = 0; d < D; d++)
        {
            m->elts[d][0] = m_t->elts[0][d] = x->elts[n][d] - mean->elts[d];
        }

        multiply_matrices_d(&M, m, m_t);
        add_matrices_d(&cov, cov, M);
    }

    multiply_matrix_by_scalar_d(&cov, cov, 1.0/(double)(N - 1));

    free_vector_d(mean);
    free_matrix_d(m);
    free_matrix_d(m_t);
    free_matrix_d(M);
}

void mv_sample_stats_f
(
    Vector_f**      mean_out, 
    Matrix_f**      cov_out,
    const Matrix_f* x
)
{
    uint32_t n, N;
    uint32_t d, D;
    Vector_f* mean;
    Matrix_f* cov;
    Matrix_f* m   = 0;
    Matrix_f* m_t = 0;
    Matrix_f* M   = 0;

    D = x->num_cols;
    N = x->num_rows;

    create_zero_vector_f(mean_out, D);
    mean = *mean_out;

    for (n = 0; n < N; n++)
    {
        for (d = 0; d < D; d++)
        {
            mean->elts[ d ] += x->elts[ n ][ d ];
        }
    }

    multiply_vector_by_scalar_f(&mean, mean, 1.0 / (float)N);

    create_zero_matrix_f(cov_out, D, D);
    cov = *cov_out;

    create_matrix_f(&m, D, 1);
    create_matrix_f(&m_t, 1, D);
    create_matrix_f(&M, D, D);

    for (n = 0; n < N; n++)
    {
        for (d = 0; d < D; d++)
        {
            m->elts[d][0] = m_t->elts[0][d] = x->elts[n][d] - mean->elts[d];
        }

        multiply_matrices_f(&M, m, m_t);
        add_matrices_f(&cov, cov, M);
    }

    if (N > 1)
    {
        multiply_matrix_by_scalar_f(&cov, cov, 1.0/(double)(N - 1));
    }

    free_matrix_f(m);
    free_matrix_f(m_t);
    free_matrix_f(M);
}

void mv_sample_stats_d
(
    Vector_d**      mean_out, 
    Matrix_d**      cov_out, 
    const Matrix_d* x
)
{
    uint32_t n, N;
    uint32_t d, D;
    Vector_d* mean;
    Matrix_d* cov;
    Matrix_d* m   = 0;
    Matrix_d* m_t = 0;
    Matrix_d* M   = 0;

    D = x->num_cols;
    N = x->num_rows;

    create_zero_vector_d(mean_out, D);
    mean = *mean_out;

    for (n = 0; n < N; n++)
    {
        for (d = 0; d < D; d++)
        {
            mean->elts[ d ] += x->elts[ n ][ d ];
        }
    }

    multiply_vector_by_scalar_d(&mean, mean, 1.0 / (double)N);

    create_zero_matrix_d(cov_out, D, D);
    cov = *cov_out;

    create_matrix_d(&m, D, 1);
    create_matrix_d(&m_t, 1, D);
    create_matrix_d(&M, D, D);

    for (n = 0; n < N; n++)
    {
        for (d = 0; d < D; d++)
        {
            m->elts[d][0] = m_t->elts[0][d] = x->elts[n][d] - mean->elts[d];
        }

        multiply_matrices_d(&M, m, m_t);
        add_matrices_d(&cov, cov, M);
    }

    if (N > 1)
    {
        multiply_matrix_by_scalar_d(&cov, cov, 1.0/(double)(N - 1));
    }

    free_matrix_d(m);
    free_matrix_d(m_t);
    free_matrix_d(M);
}