estimate.c

#include <math.h>
#include <stdlib.h>
#include "tabledesign.h"

/**
 * Computes the autocorrelation of a vector. More precisely, it computes the
 * dot products of vec[i:] and vec[:-i] for i in [0, k). Unused.
 *
 * See https://en.wikipedia.org/wiki/Autocorrelation.
 */
void acf(double *vec, int n, double *out, int k)
{
    int i, j;
    double sum;
    for (i = 0; i < k; i++)
    {
        sum = 0.0;
        for (j = 0; j < n - i; j++)
        {
            sum += vec[j + i] * vec[j];
        }
        out[i] = sum;
    }
}

// https://en.wikipedia.org/wiki/Durbin%E2%80%93Watson_statistic ?
// "detects the presence of autocorrelation at lag 1 in the residuals (prediction errors)"
int durbin(double *arg0, int n, double *arg2, double *arg3, double *outSomething)
{
    int i, j;
    double sum, div;
    int ret;

    arg3[0] = 1.0;
    div = arg0[0];
    ret = 0;

    for (i = 1; i <= n; i++)
    {
        sum = 0.0;
        for (j = 1; j <= i-1; j++)
        {
            sum += arg3[j] * arg0[i - j];
        }

        arg3[i] = (div > 0.0 ? -(arg0[i] + sum) / div : 0.0);
        arg2[i] = arg3[i];

        if (fabs(arg2[i]) > 1.0)
        {
            ret++;
        }

        for (j = 1; j < i; j++)
        {
            arg3[j] += arg3[i - j] * arg3[i];
        }

        div *= 1.0 - arg3[i] * arg3[i];
    }
    *outSomething = div;
    return ret;
}

void afromk(double *in, double *out, int n)
{
    int i, j;
    out[0] = 1.0;
    for (i = 1; i <= n; i++)
    {
        out[i] = in[i];
        for (j = 1; j <= i - 1; j++)
        {
            out[j] += out[i - j] * out[i];
        }
    }
}

int kfroma(double *in, double *out, int n)
{
    int i, j;
    double div;
    double temp;
    double *next;
    int ret;

    ret = 0;
    next = malloc((n + 1) * sizeof(double));

    out[n] = in[n];
    for (i = n - 1; i >= 1; i--)
    {
        for (j = 0; j <= i; j++)
        {
            temp = out[i + 1];
            div = 1.0 - (temp * temp);
            if (div == 0.0)
            {
                free(next);
                return 1;
            }
            next[j] = (in[j] - in[i + 1 - j] * temp) / div;
        }

        for (j = 0; j <= i; j++)
        {
            in[j] = next[j];
        }

        out[i] = next[i];
        if (fabs(out[i]) > 1.0)
        {
            ret++;
        }
    }

    free(next);
    return ret;
}

void rfroma(double *arg0, int n, double *arg2)
{
    int i, j;
    double **mat;
    double div;

    mat = malloc((n + 1) * sizeof(double*));
    mat[n] = malloc((n + 1) * sizeof(double));
    mat[n][0] = 1.0;
    for (i = 1; i <= n; i++)
    {
        mat[n][i] = -arg0[i];
    }

    for (i = n; i >= 1; i--)
    {
        mat[i - 1] = malloc(i * sizeof(double));
        div = 1.0 - mat[i][i] * mat[i][i];
        for (j = 1; j <= i - 1; j++)
        {
            mat[i - 1][j] = (mat[i][i - j] * mat[i][i] + mat[i][j]) / div;
        }
    }

    arg2[0] = 1.0;
    for (i = 1; i <= n; i++)
    {
        arg2[i] = 0.0;
        for (j = 1; j <= i; j++)
        {
            arg2[i] += mat[i][j] * arg2[i - j];
        }
    }

    free(mat[n]);
    for (i = n; i > 0; i--)
    {
        free(mat[i - 1]);
    }
    free(mat);
}

double model_dist(double *arg0, double *arg1, int n)
{
    double *sp3C;
    double *sp38;
    double ret;
    int i, j;

    sp3C = malloc((n + 1) * sizeof(double));
    sp38 = malloc((n + 1) * sizeof(double));
    rfroma(arg1, n, sp3C);

    for (i = 0; i <= n; i++)
    {
        sp38[i] = 0.0;
        for (j = 0; j <= n - i; j++)
        {
            sp38[i] += arg0[j] * arg0[i + j];
        }
    }

    ret = sp38[0] * sp3C[0];
    for (i = 1; i <= n; i++)
    {
        ret += 2 * sp3C[i] * sp38[i];
    }

    free(sp3C);
    free(sp38);
    return ret;
}

// compute autocorrelation matrix?
void acmat(short *in, int n, int m, double **out)
{
    int i, j, k;
    for (i = 1; i <= n; i++)
    {
        for (j = 1; j <= n; j++)
        {
            out[i][j] = 0.0;
            for (k = 0; k < m; k++)
            {
                out[i][j] += in[k - i] * in[k - j];
            }
        }
    }
}

// compute autocorrelation vector?
void acvect(short *in, int n, int m, double *out)
{
    int i, j;
    for (i = 0; i <= n; i++)
    {
        out[i] = 0.0;
        for (j = 0; j < m; j++)
        {
            out[i] -= in[j - i] * in[j];
        }
    }
}

/**
 * Replaces a real n-by-n matrix "a" with the LU decomposition of a row-wise
 * permutation of itself.
 *
 * Input parameters:
 * a: The matrix which is operated on. 1-indexed; it should be of size
 *    (n+1) x (n+1), and row/column index 0 is not used.
 * n: The size of the matrix.
 *
 * Output parameters:
 * indx: The row permutation performed. 1-indexed; it should be of size n+1,
 *       and index 0 is not used.
 * d: the determinant of the permutation matrix.
 *
 * Returns 1 to indicate failure if the matrix is singular or has zeroes on the
 * diagonal, 0 on success.
 *
 * Derived from ludcmp in "Numerical Recipes in C: The Art of Scientific Computing",
 * with modified error handling.
 */
int lud(double **a, int n, int *indx, int *d)
{
    int i,imax,j,k;
    double big,dum,sum,temp;
    double min,max;
    double *vv;

    vv = malloc((n + 1) * sizeof(double));
    *d=1;
    for (i=1;i<=n;i++) {
        big=0.0;
        for (j=1;j<=n;j++)
            if ((temp=fabs(a[i][j])) > big) big=temp;
        if (big == 0.0) return 1;
        vv[i]=1.0/big;
    }
    for (j=1;j<=n;j++) {
        for (i=1;i<j;i++) {
            sum=a[i][j];
            for (k=1;k<i;k++) sum -= a[i][k]*a[k][j];
            a[i][j]=sum;
        }
        big=0.0;
        for (i=j;i<=n;i++) {
            sum=a[i][j];
            for (k=1;k<j;k++)
                sum -= a[i][k]*a[k][j];
            a[i][j]=sum;
            if ( (dum=vv[i]*fabs(sum)) >= big) {
                big=dum;
                imax=i;
            }
        }
        if (j != imax) {
            for (k=1;k<=n;k++) {
                dum=a[imax][k];
                a[imax][k]=a[j][k];
                a[j][k]=dum;
            }
            *d = -(*d);
            vv[imax]=vv[j];
        }
        indx[j]=imax;
        if (a[j][j] == 0.0) return 1;
        if (j != n) {
            dum=1.0/(a[j][j]);
            for (i=j+1;i<=n;i++) a[i][j] *= dum;
        }
    }
    free(vv);

    min = 1e10;
    max = 0.0;
    for (i = 1; i <= n; i++)
    {
        temp = fabs(a[i][i]);
        if (temp < min) min = temp;
        if (temp > max) max = temp;
    }
    return min / max < 1e-10 ? 1 : 0;
}

/**
 * Solves the set of n linear equations Ax = b, using LU decomposition
 * back-substitution.
 *
 * Input parameters:
 * a: The LU decomposition of a matrix, created by "lud".
 * n: The size of the matrix.
 * indx: Row permutation vector, created by "lud".
 * b: The vector b in the equation. 1-indexed; is should be of size n+1, and
 *    index 0 is not used.
 *
 * Output parameters:
 * b: The output vector x. 1-indexed.
 *
 * From "Numerical Recipes in C: The Art of Scientific Computing".
 */
void lubksb(double **a, int n, int *indx, double *b)
{
    int i,ii=0,ip,j;
    double sum;

    for (i=1;i<=n;i++) {
        ip=indx[i];
        sum=b[ip];
        b[ip]=b[i];
        if (ii)
            for (j=ii;j<=i-1;j++) sum -= a[i][j]*b[j];
        else if (sum) ii=i;
        b[i]=sum;
    }
    for (i=n;i>=1;i--) {
        sum=b[i];
        for (j=i+1;j<=n;j++) sum -= a[i][j]*b[j];
        b[i]=sum/a[i][i];
    }
}