package jhpro.stat;
import java.util.*;
import org.apache.commons.math3.util.FastMath;
/**
* Common statistical methods.
*
* @author S.Chekanov.
*
*/
public class StatUtils {
private static Random _random;
/**
* This static class does the autocorrelation for the given float vector. It
* uses the function: autocorrelation[lag] = SUM(n=0,N)
* (signal[n]*signal[n+lag]) with n is the input sample index, and
* 0signal at different
* times n are the same values of the same signal delays
* by lag samples determines the magnitude of
* autocorrelation[lag]. The output of an autocorrelation shows the
* magnitude for different lag times. Note that the array length of the
* signal vector must have a length from at least window length + window
* shift (signal.length >= N+lag-max). If the length is smaller a null float
* array will be given back.
*
* @param float[] signal - the signal vector.
* @param int windowlength - the length of correlation build with the
* summaration loop. Corresponds to lag-max
* @param int windowshift - the length of the summaration lopp. Corresponds
* to N
* @return float[] - the autocorrelated signal vector.
*/
public static double[] autoCorrelate(double[] signal, int windowlength,
int windowshift) {
int signallenth = signal.length;
if (signallenth < windowlength + windowshift)
return null;
double[] autocorrelation = new double[windowlength];
// loop over the magnitude of the autocorrelation
for (int lag = 0; lag < windowlength; lag++) {
// loop over the sum
for (int n = 0; n < windowshift; n++) {
autocorrelation[lag] += signal[n] * signal[n + lag];
}
}
return autocorrelation;
}
/**
* Compute average of all values in vector v.
*
* @param v
* a vector of float[]
* @return the average of all elements in v.
*/
public static double avg(double[] v) {
return (float) (sum(v) / v.length);
}
/**
* sum(1) computes summation of all values in vector v.
*
* @param v
* a vector of float[]
* @return the average of all elements in v.
*/
public static double sum(float[] v) {
double aggr = 0;
for (int j = 0; j < v.length; j++)
aggr += v[j];
return aggr;
}
/**
* sum(1) computes summation of all values in matrix m.
*
* @param v
* a vector of float[]
* @return the average of all elements in v.
*/
public static double sum(float[][] m) {
double aggr = 0;
for (int j = 0; j < m.length; j++)
for (int i = 0; i < m[0].length; i++)
aggr += m[j][i];
return aggr;
}
/**
* sum(2) computes the sum of the elements of vector v only where the mask
* value is not 0.
*
* @param v
* a float[] vector.
* @param mask
* a float[] with a value of ! 0 for all elements of v that are
* valid.
* @return the avg of the masked elements in v.
*/
public static double sum(float[] v, float[] mask) {
double agr = 0;
for (int j = 0; j < v.length; j++) {
if (mask[j] != 0)
agr += v[j];
}
return agr;
}
/**
* Compute average of all values in vector v.
*
* @param v
* a vector of float[]
* @return the average of all elements in v.
*/
public static double sum(double[] v) {
double aggr = 0;
for (int j = 0; j < v.length; j++)
aggr += v[j];
return aggr;
}
/**
* Compute stdev (SQRT(var)) of all values in vector v.
*
* @param v
* a vector of float[]
*/
public static double stdev(double[] v) {
float stdev = (float) FastMath.sqrt(var(v));
return stdev;
}
/**
* Compute the correlation of a vector with another vector b. The mean of a
* and b is subtracted before further processing. a and b therefore do not
* have to be zero-mean.
*
* @param a
* a float[] vector.
* @param b
* a float[] vector of same length.
* @return r, the correlation coefficient.
*/
public static double correl(double[] a, double[] b) {
int iN = FastMath.min(a.length, b.length);
double avga = avg(a);
double avgb = avg(b);
double agr = 0;
double vara = 0;
double varb = 0;
for (int i = 0; i < iN; i++) {
agr += (a[i] - avga) * (b[i] - avgb);
vara += FastMath.pow((a[i] - avga), 2);
varb += FastMath.pow((b[i] - avgb), 2);
}
double r;
if (vara == 0 || varb == 0)
r = 0;
else
r = (agr / FastMath.sqrt(vara * varb));
if (Double.isNaN(r))
System.out.println("NAN->" + agr + " " + vara + " " + varb);
return (float) r;
}
/**
* Compute variance of all elements in vector v.
*
* @param v
* a vector of float[]
*/
public static double var(double[] v) {
double avg = avg(v);
double var = 0;
double eps = 0;
int n = 0;
for (int j = 0; j < v.length; j++) {
double d = v[j] - avg;
eps += d;
var += d * d;
n++;
}
var = (var - eps * eps / n) / (n - 1);
return var;
}
/**
* Compute residuals of b versus a.
*
* @param a
* a float[] vector
* @param b
* a float[] vector.
* @return a float[] that contains the residuals for each b-a.
*/
public static float[] residuals(float[] a, float[] b)
throws IllegalArgumentException {
float[] res = new float[FastMath.min(a.length, b.length)];
for (int i = 0; i < res.length; i++)
res[i] = b[i] - a[i];
return res;
}
/**
* Computes autocovariance up to maxk for the time series data.
*/
public static double[] autocov(double[] data, int maxk) {
double[] ac = new double[maxk];
data = zeromean(data);
data = linearcorrect(data);
double mu = avg(data);
for (int k = 0; k < maxk; k++) {
for (int t = 0; t < data.length - k; t++) {
ac[k] += (data[t] - mu) * (data[t + k] - mu);
}
ac[k] /= data.length - k;
}
return ac;
}
/**
* Corrects data for a linear trend.
*/
protected static double[] linearcorrect(double[] data) {
double sumx = 0;
double sumy = 0;
double sumxy = 0;
double sumxx = 0;
double sumyy = 0;
for (int k = 0; k < data.length; k++) {
sumx += k;
sumy += data[k];
sumxx += k * k;
sumyy += data[k] * data[k];
sumxy += data[k] * k;
}
double Sxx = sumxx - sumx * sumx / data.length;
double yy = sumyy - sumy * sumy / data.length;
double Sxy = sumxy - sumx * sumy / data.length;
double b = Sxy / Sxx;
double a = (sumy - b * sumx) / data.length;
for (int k = 0; k < data.length; k++) {
data[k] = data[k] - b * k;
}
return data;
}
/**
* Computes autocorrelation up to maxk for the time series data.
*/
public static double[] autocorr(double[] data, int maxk) {
double[] acov = autocov(data, maxk);
double[] acorr = new double[acov.length];
for (int k = 0; k < acov.length; k++) {
acorr[k] = acov[k] / acov[0];
}
return acorr;
}
/**
* Zero the mean of all elements of v.
*
* @param v
* a vector.
* @return a vector with zero mean and same variance as v.
*/
public static double[] zeromean(double[] v) {
double[] n = new double[v.length];
double avg = sum(v) / v.length;
// Make zero mean.
for (int j = 0; j < v.length; j++)
n[j] = v[j] - avg;
return n;
}
/**
* Return a new array that is the cross-correlation of the two argument
* arrays, starting and ending at user-specified lag values. The output
* array will have length (endLag - startLag + 1). The first element of the
* output will have the cross-correlation at a lag of startLag. The last
* element of the output will have the cross-correlation at a lag of endLag.
*
* @param x
* The first array of doubles.
* @param y
* The second array of doubles.
* @param N
* An integer indicating the number of samples to sum over.
* @param startLag
* An int indicating at which lag to start (may be negative).
* @param endLag
* An int indicating at which lag to end.
* @return A new array of doubles.
*/
public static final double[] crossCorrelation(double[] x, double[] y,
int N, int startLag, int endLag) {
int outputLength = endLag - startLag + 1;
double[] returnValue = new double[outputLength];
for (int lag = startLag; lag <= endLag; lag++) {
// Find the most efficient and correct place to start the summation
int start = FastMath.max(0, -lag);
// Find the most efficient and correct place to end the summation
int limit = FastMath.min(x.length, N);
limit = FastMath.min(limit, y.length - lag);
double sum = 0.0;
for (int i = start; i < limit; i++) {
sum += (x[i] * y[i + lag]);
}
returnValue[lag - startLag] = sum;
}
return returnValue;
}
/**
* Return the cross-correlation of two arrays at a certain lag value. The
* cross-correlation is defined by : Rxy[d] = sum of i = 0 to N - 1 of x[i]
* * y[i + d]
*
* @param x
* The first array of doubles.
* @param y
* The second array of doubles.
* @param N
* An integer indicating the number of samples to sum over. This
* must be non-negative, but large numbers are ok because this
* routine will not overrun reading of the arrays.
* @param lag
* An integer indicating the lag value (may be negative).
* @return A double, Rxy[lag].
*/
public static double crossCorrelationAt(double[] x, double[] y, int N,
int lag) {
// Find the most efficient and correct place to start the summation
int start = FastMath.max(0, -lag);
// Find the most efficient and correct place to end the summation
int limit = FastMath.min(x.length, N);
limit = FastMath.min(limit, y.length - lag);
double sum = 0.0;
for (int i = start; i < limit; i++) {
sum += (x[i] * y[i + lag]);
}
return sum;
}
/**
* Return a new array of Bernoulli random variables with a given probability
* of success p. On success, the random variable has value 1.0; on failure
* the random variable has value 0.0.
*
* @param p
* The probability, which should be a double between 0.0 and 1.0.
* The probability is compared to the output of
* java.lang.Random.nextDouble(). If the output is less than p,
* then the array element will be 1.0. If the output is greater
* than or equal to p, then the array element will be 0.0.
* @param N
* The number of elements to allocate.
* @return The array of Bernoulli random variables.
*/
public static final double[] randomBernoulli(double p, int N) {
double[] returnValue = new double[N];
if (_random == null) {
_random = new Random();
}
for (int i = 0; i < N; i++) {
returnValue[i] = (_random.nextDouble() < p) ? 1.0 : 0.0;
}
return returnValue;
}
/**
* Return a new array of exponentially distributed doubles with parameter
* lambda. The number of elements to allocate is given by N.
*
* @param lambda
* The lambda, which may not by 0.0.
* @param N
* The number of elements to allocate.
* @return The array of exponential random variables.
*/
public static final double[] randomExponential(double lambda, int N) {
double[] returnValue = new double[N];
if (_random == null) {
_random = new Random();
}
for (int i = 0; i < N; i++) {
double r;
do {
r = _random.nextDouble();
} while (r == 0.0);
returnValue[i] = -FastMath.log(r) / lambda;
}
return returnValue;
}
/**
* Return a new array of Gaussian distributed doubles with a given mean and
* standard deviation. The number of elements to allocate is given by N.
* This algorithm is from [1].
*
* @param mean
* The mean of the new array.
* @param standardDeviation
* The standard deviation of the new array.
* @param N
* The number of elements to allocate.
* @return The array of random Gaussian variables.
*/
public static final double[] randomGaussian(double mean,
double standardDeviation, int N) {
double[] returnValue = new double[N];
if (_random == null) {
_random = new Random();
}
for (int i = 0; i < N; i++) {
returnValue[i] = mean
+ (_random.nextGaussian() * standardDeviation);
}
return returnValue;
}
/**
* Return a new array of Poisson random variables (as doubles) with a given
* mean. The number of elements to allocate is given by N. This algorithm is
* from [1].
*
* @param mean
* The mean of the new array.
* @param N
* The number of elements to allocate.
* @return The array of random Poisson variables.
*/
public static final double[] randomPoisson(double mean, int N) {
double[] returnValue = new double[N];
if (_random == null) {
_random = new Random();
}
for (int i = 0; i < N; i++) {
double j;
double u;
double p;
double f;
j = 0.0;
f = p = FastMath.exp(-mean);
u = _random.nextDouble();
while (f <= u) {
p *= (mean / (j + 1.0));
f += p;
j += 1.0;
}
returnValue[i] = j;
}
return returnValue;
}
/**
* Return a new array of uniformly distributed doubles ranging from a to b.
* The number of elements to allocate is given by N.
*
* @param a
* A double indicating the lower bound.
* @param b
* A double indicating the upper bound.
* @param N
* An int indicating how many elements to generate.
* @return A new array of doubles.
*/
public static double[] randomUniform(double a, double b, int N) {
double range = b - a;
double[] returnValue = new double[N];
if (_random == null) {
_random = new Random();
}
for (int i = 0; i < N; i++) {
returnValue[i] = (_random.nextDouble() * range) + a;
}
return returnValue;
}
}
Ads help maintain this website.