Example usage for org.apache.commons.math3.special Gamma regularizedGammaQ

List of usage examples for org.apache.commons.math3.special Gamma regularizedGammaQ

Introduction

In this page you can find the example usage for org.apache.commons.math3.special Gamma regularizedGammaQ.

Prototype

public static double regularizedGammaQ(double a, double x) 

Source Link

Document

Returns the regularized gamma function Q(a, x) = 1 - P(a, x).

Usage

From source file:com.intel.diceros.test.securerandom.DRNGTest.java

/**
 * Frequency Test with a Block//  www .  ja v  a2s .com
 * <p>
 * The focus of the test is the proportion of ones within M-bit
 * blocks. The purpose of this test is to determine whether the
 * frequency of ones in an M-bit block is approximately block_size/2,
 * as would be expected under an assumption of randomness. For block
 * size = 1, this test degenerates to test 1, frequency test for RNG.
 */
private void testBlockFrequency(final int[] epsilon, final int M) {
    final int n = epsilon.length;
    final int N = n / M;

    double sum = 0.0;

    for (int i = 0; i < N; i++) {
        int blockSum = 0;
        for (int j = 0; j < M; j++) {
            blockSum += epsilon[j + i * M];
        }
        double pi = (double) blockSum / (double) M;
        double v = pi - 0.5;
        sum += v * v;
    }
    double chi_squared = 4.0 * M * sum;
    double p_value = Gamma.regularizedGammaQ(N / 2, chi_squared / 2);

    assertTrue("RNG test failed, test block frequency.", p_value >= 0.01);
}

From source file:com.intel.diceros.test.securerandom.DRNGTest.java

/**
 * Test for the Longest Run of Ones in a Block
 * <p>//from  w  ww  .  ja v  a  2 s.c  o  m
 * The focus of the test is the longest run of ones within M-bit blocks.
 * The purpose of this test is to determine whether the length of the
 * longest run of ones within the tested sequence is consistent with the
 * length of the longest run of ones that would be expected in a random
 * sequence. Note that an irregularity in the exported length of the
 * longest run of ones implies that there is also an irregularity in the
 * expected length of the longest run of zeros. Therefor, ony a test for
 * ones is necessary.
 */
private void testLongestRunOfOnes(final int[] epsilon) {
    final int n = epsilon.length;

    int k, m, v_n_obs, run;
    int[] v = new int[7];
    int[] nu = { 0, 0, 0, 0, 0, 0, 0 };
    double[] pi = new double[7];

    if (n < 128)
        return;

    if (n < 6272) {
        k = 3;
        m = 8;
        v[0] = 1;
        v[1] = 2;
        v[2] = 3;
        v[3] = 4;
        pi[0] = 0.21484375;
        pi[1] = 0.3671875;
        pi[2] = 0.23046875;
        pi[3] = 0.1875;
    } else if (n < 750000) {
        k = 5;
        m = 128;
        v[0] = 4;
        v[1] = 5;
        v[2] = 6;
        v[3] = 7;
        v[4] = 8;
        v[5] = 9;
        pi[0] = 0.1174035788;
        pi[1] = 0.242955959;
        pi[2] = 0.249363483;
        pi[3] = 0.17517706;
        pi[4] = 0.102701071;
        pi[5] = 0.112398847;
    } else {
        k = 6;
        m = 10000;
        v[0] = 10;
        v[1] = 11;
        v[2] = 12;
        v[3] = 13;
        v[4] = 14;
        v[5] = 15;
        v[6] = 16;
        pi[0] = 0.0882;
        pi[1] = 0.2092;
        pi[2] = 0.2483;
        pi[3] = 0.1933;
        pi[4] = 0.1208;
        pi[5] = 0.0675;
        pi[6] = 0.0727;
    }

    int N = n / m;
    for (int i = 0; i < N; i++) {
        v_n_obs = 0;
        run = 0;
        for (int j = 0; j < m; j++) {
            if (epsilon[i * m + j] == 1) {
                if (++run > v_n_obs)
                    v_n_obs = run;
            } else
                run = 0;
        }
        if (v_n_obs < v[0])
            nu[0]++;
        for (int j = 0; j <= k; j++) {
            if (v_n_obs == v[j])
                nu[j]++;
        }
        if (v_n_obs > v[k])
            nu[k]++;
    }

    double chi2 = 0.0;
    for (int i = 0; i < k; i++)
        chi2 += ((nu[i] - N * pi[i]) * (nu[i] - N * pi[i])) / (N * pi[i]);

    double p_value = Gamma.regularizedGammaQ(k / 2.0, chi2 / 2.0);
    assertTrue("RNG test failed, test longest run of ones.", p_value >= 0.01);
}

From source file:com.intel.diceros.test.securerandom.DRNGTest.java

/**
 * Test Non-overlapping Template Matching Test
 * <p>//from  ww  w  . j  av  a  2s.  co m
 * The focus of this test is the number of occurrences of pre-specified target
 * strings. The purpose of this test is to detect generators that produce too
 * many occurrences of a given non-periodic (aperiodic) pattern. For this test
 * an m-bit window is used to search for a specific m-bit pattern. If the
 * pattern is not found, the window slides one bit position. If the pattern is
 * found, the window is reset to the bit after the found pattern, and the search
 * resumes.
 */
private void testNonOverlappingTemplateMatching(final int[] epsilon, final int m) throws IOException {

    int[] numOfTemplates = { 0, 0, 2, 4, 6, 12, 20, 40, 74, 148, 284, 568, 1116, 2232, 4424, 8848, 17622, 35244,
            70340, 140680, 281076, 562152 };
    final int maxNumOfTemplates = numOfTemplates.length;
    numOfTemplates = Arrays.copyOf(numOfTemplates, 100);

    int i, j, jj, k, match, skip, M, N, K = 5, n = epsilon.length;
    N = 8;
    M = n / N;
    int w_obs;
    int[] nu = new int[6], wj = new int[N], sequence = new int[m];
    double sum, chi2, p_value, lambda, varWj;
    double[] pi = new double[6];

    lambda = (M - m + 1) / Math.pow(2, m);
    assertTrue("Lambda not being positive!", lambda > 0);

    varWj = M * (1.0 / Math.pow(2.0, m) - (2.0 * m - 1.0) / Math.pow(2.0, 2.0 * m));

    if (numOfTemplates[m] < maxNumOfTemplates)
        skip = 1;
    else
        skip = numOfTemplates[m] / maxNumOfTemplates;
    numOfTemplates[m] /= skip;

    sum = 0.0;
    for (i = 0; i < 2; i++) {
        pi[i] = Math.exp(-lambda + i * Math.log(lambda) - Gamma.logGamma(i + 1));
        sum += pi[i];
    }
    pi[0] = sum;
    for (i = 2; i < K; i++) {
        pi[i - 1] = Math.exp(-lambda + i * Math.log(lambda) - Gamma.logGamma(i + 1));
        sum += pi[i - 1];
    }
    pi[K] = 1 - sum;

    final String templateName = "template" + m;
    Reader reader = null;
    try {
        reader = getReaders(templateName).get(0);
        for (jj = 0; jj < Math.min(maxNumOfTemplates, numOfTemplates[m]); jj++) {
            sum = 0;
            for (k = 0; k < m; k++) {
                int b = reader.read();
                if (b == -1)
                    assertTrue("Template issue.", false);
                sequence[k] = b - 48;
            }
            for (k = 0; k <= K; k++)
                nu[k] = 0;
            for (i = 0; i < N; i++) {
                w_obs = 0;
                for (j = 0; j < M - m + 1; j++) {
                    match = 1;
                    for (k = 0; k < m; k++) {
                        if (sequence[k] != epsilon[i * M + j + k]) {
                            match = 0;
                            break;
                        }
                    }
                    if (match == 1)
                        w_obs++;
                }
                wj[i] = w_obs;
            }
            //            sum = 0;
            chi2 = 0.0; /* Compute Chi Square */
            for (i = 0; i < N; i++) {
                chi2 += Math.pow(((double) wj[i] - lambda) / Math.pow(varWj, 0.5), 2);
            }
            p_value = Gamma.regularizedGammaQ(N / 2.0, chi2 / 2.0);
            assertTrue("RNG test failed, test non-overlapping template matching.", p_value >= 0.01);
            if (skip > 1)
                for (i = 0; i < (skip - 1) * 2 * m; i++)
                    reader.read();
        }
    } finally {
        if (reader != null)
            reader.close();
    }
}

From source file:com.intel.diceros.test.securerandom.DRNGTest.java

/**
 * Overlapping Template Matching Test//from   w  ww . j a v a  2 s .c  om
 * <p>
 * The focus of the Overlapping Template Matching test is the number of
 * occurrences of pre-specified target strings. Both this test uses an
 * m-bit window to search for a specific m-bit pattern. If the pattern
 * is not found, the window slides one bit position.
 */
private void testOverlappingTemplateMatchings(final int[] epsilon, final int m) {
    int i, k, match;
    double w_obs, eta, sum, chi2, p_value, lambda;
    int M, N, j, K = 5;
    int[] nu = { 0, 0, 0, 0, 0, 0 }, sequence = new int[m];
    double[] pi = { 0.143783, 0.139430, 0.137319, 0.124314, 0.106209, 0.348945 };
    final int n = epsilon.length;
    M = 1032;
    N = n / M;

    for (i = 0; i < m; i++)
        sequence[i] = 1;

    lambda = (double) (M - m + 1) / Math.pow(2, m);
    eta = lambda / 2.0;
    sum = 0.0;
    for (i = 0; i < K; i++) { /* Compute Probabilities */
        pi[i] = pr(i, eta);
        sum += pi[i];
    }
    pi[K] = 1 - sum;

    for (i = 0; i < N; i++) {
        w_obs = 0;
        for (j = 0; j < M - m + 1; j++) {
            match = 1;
            for (k = 0; k < m; k++) {
                if (sequence[k] != epsilon[i * M + j + k])
                    match = 0;
            }
            if (match == 1)
                w_obs++;
        }
        if (w_obs <= 4)
            nu[(int) w_obs]++;
        else
            nu[K]++;
    }
    sum = 0;
    chi2 = 0.0; /* Compute Chi Square */
    for (i = 0; i < K + 1; i++) {
        chi2 += Math.pow((double) nu[i] - (double) N * pi[i], 2) / ((double) N * pi[i]);
        sum += nu[i];
    }
    p_value = Gamma.regularizedGammaQ(K / 2.0, chi2 / 2.0);
    assertTrue("RNG test failed, test overlapping template matching.", p_value >= 0.01);
}

From source file:com.intel.diceros.test.securerandom.DRNGTest.java

/**
 * Linear Complexity Test/*www  .j a  va  2s.c  o m*/
 * <p>
 * The focus of this test is the length of a linear feedback
 * shift register (LFSR). The purpose of this test is to determine
 * whether or not the sequence is complex enough to be considered
 * random. Random sequences are characterized by longer LFSRs. An
 * LFSR that is too short implies non-randomness.
 */
private void testLinearComplexity(final int[] epsilon, int M) {
    int i, ii, j, d, N, L, m, N_, parity, sign, K = 6, n = epsilon.length;
    double p_value, T_, mean, chi2;
    double[] pi = { 0.01047, 0.03125, 0.12500, 0.50000, 0.25000, 0.06250, 0.020833 }, nu = new double[7];
    int[] T = new int[M], P = new int[M], B_ = new int[M], C = new int[M];

    N = (int) Math.floor(n / M);
    for (i = 0; i < K + 1; i++)
        nu[i] = 0.00;
    for (ii = 0; ii < N; ii++) {
        for (i = 0; i < M; i++) {
            B_[i] = 0;
            C[i] = 0;
            T[i] = 0;
            P[i] = 0;
        }
        L = 0;
        m = -1;
        d = 0;
        C[0] = 1;
        B_[0] = 1;

        /* DETERMINE LINEAR COMPLEXITY */
        N_ = 0;
        while (N_ < M) {
            d = epsilon[ii * M + N_];
            for (i = 1; i <= L; i++)
                d += C[i] * epsilon[ii * M + N_ - i];
            d = d % 2;
            if (d == 1) {
                for (i = 0; i < M; i++) {
                    T[i] = C[i];
                    P[i] = 0;
                }
                for (j = 0; j < M; j++)
                    if (B_[j] == 1)
                        P[j + N_ - m] = 1;
                for (i = 0; i < M; i++)
                    C[i] = (C[i] + P[i]) % 2;
                if (L <= N_ / 2) {
                    L = N_ + 1 - L;
                    m = N_;
                    for (i = 0; i < M; i++)
                        B_[i] = T[i];
                }
            }
            N_++;
        }
        if ((parity = (M + 1) % 2) == 0)
            sign = -1;
        else
            sign = 1;
        mean = M / 2.0 + (9.0 + sign) / 36.0 - 1.0 / Math.pow(2, M) * (M / 3.0 + 2.0 / 9.0);
        if ((parity = M % 2) == 0)
            sign = 1;
        else
            sign = -1;
        T_ = sign * (L - mean) + 2.0 / 9.0;

        if (T_ <= -2.5)
            nu[0]++;
        else if (T_ > -2.5 && T_ <= -1.5)
            nu[1]++;
        else if (T_ > -1.5 && T_ <= -0.5)
            nu[2]++;
        else if (T_ > -0.5 && T_ <= 0.5)
            nu[3]++;
        else if (T_ > 0.5 && T_ <= 1.5)
            nu[4]++;
        else if (T_ > 1.5 && T_ <= 2.5)
            nu[5]++;
        else
            nu[6]++;
    }
    chi2 = 0.00;
    for (i = 0; i < K + 1; i++)
        chi2 += Math.pow(nu[i] - N * pi[i], 2) / (N * pi[i]);
    p_value = Gamma.regularizedGammaQ(K / 2.0, chi2 / 2.0);

    assertTrue("RNG test failed, test linear complexity.", p_value >= 0.01);
}

From source file:com.intel.diceros.test.securerandom.DRNGTest.java

/**
 * Serial Test// w ww. java2 s . com
 * <p>
 * The focus of this test is the frequency of all possible overlapping m-bit
 * patterns across the entire sequence. The purpose of this test is to determine
 * whether the number of occurrences of the 2mm-bit overlapping patterns is
 * approximately the same as would be expected for a random sequence. Random
 * sequences have uniformity; that is, every m-bit pattern has the same chance
 * of appearing as every other m-bit pattern. Note that for m = 1, the Serial
 * test is equivalent to the Frequency test of Section 2.1.
 */
private void testSerial(final int[] epsilon, int m) {
    double p_value1, p_value2, psim0, psim1, psim2, del1, del2;
    psim0 = psi2(epsilon, m);
    psim1 = psi2(epsilon, m - 1);
    psim2 = psi2(epsilon, m - 2);
    del1 = psim0 - psim1;
    del2 = psim0 - 2.0 * psim1 + psim2;
    p_value1 = Gamma.regularizedGammaQ(Math.pow(2, m - 1) / 2, del1 / 2.0);
    p_value2 = Gamma.regularizedGammaQ(Math.pow(2, m - 2) / 2, del2 / 2.0);

    assertTrue("RNG test failed(p_value1), test linear complexity.", p_value1 >= 0.01);
    assertTrue("RNG test failed(p_value2), test linear complexity.", p_value2 >= 0.01);
}

From source file:com.intel.diceros.test.securerandom.DRNGTest.java

/**
 * Approximate Entropy Test//from w ww  . j  a v  a2s .  c o  m
 * <p>
 * As with the Serial test of Section 2.11, the focus of this test is the frequency
 * of all possible overlapping m-bit patterns across the entire sequence.
 * The purpose of the test is to compare the frequency of overlapping blocks of
 * two consecutive/adjacent lengths (m and m+1) against the expected result for
 * a random sequence.
 */
private void testApproximateEntropy(final int[] epsilon, int m) {
    final int n = epsilon.length;
    int i, j, k, r, blockSize, seqLength, powLen, index;
    double sum, numOfBlocks, apen, chi_squared, p_value;
    double[] ApEn = new double[2];
    int[] P;

    seqLength = n;
    r = 0;
    for (blockSize = m; blockSize <= m + 1; blockSize++) {
        if (blockSize == 0) {
            ApEn[0] = 0.00;
            r++;
        } else {
            numOfBlocks = (double) seqLength;
            powLen = (int) Math.pow(2, blockSize + 1) - 1;
            P = new int[powLen];
            for (i = 1; i < powLen - 1; i++)
                P[i] = 0;
            for (i = 0; i < numOfBlocks; i++) { /* COMPUTE FREQUENCY */
                k = 1;
                for (j = 0; j < blockSize; j++) {
                    k <<= 1;
                    if ((int) epsilon[(i + j) % seqLength] == 1)
                        k++;
                }
                P[k - 1]++;
            }
            /* DISPLAY FREQUENCY */
            sum = 0.0;
            index = (int) Math.pow(2, blockSize) - 1;
            for (i = 0; i < (int) Math.pow(2, blockSize); i++) {
                if (P[index] > 0)
                    sum += P[index] * Math.log(P[index] / numOfBlocks);
                index++;
            }
            sum /= numOfBlocks;
            ApEn[r] = sum;
            r++;
        }
    }
    apen = ApEn[0] - ApEn[1];

    chi_squared = 2.0 * seqLength * (Math.log(2) - apen);
    p_value = Gamma.regularizedGammaQ(Math.pow(2, m - 1), chi_squared / 2.0);

    assertTrue("RNG test failed, test approximate entropy.", p_value >= 0.01);
}

From source file:com.intel.diceros.test.securerandom.DRNGTest.java

/**
 * Random Excursions Test//from   ww  w . java  2s .c o m
 * <p>
 * The focus of this test is the number of cycles having exactly K visits
 * in a cumulative sum random walk. The cumulative sum random walk is
 * derived from partial sums after the (0,1) sequence is transferred to
 * the appropriate (-1, +1) sequence. A cycle of a random walk consists
 * of a sequence of steps of unit length taken at random that begin at
 * and return to the origin. The purpose of this test is to determine if
 * the number of visits to a particular state within a cycle deviates
 * from what one would expect for a random sequence. This test is actually
 * a series of eight tests (and conclusions), one test and conclusion
 * for each of the states: -4, -3, -2, -1 and +1, +2, +3, +4.
 */
public void testRandomExcursions(final int[] epsilon) {
    final int n = epsilon.length;
    int[] stateX = { -4, -3, -2, -1, 1, 2, 3, 4 };
    int[] S_k = new int[n];
    int cycleMaxLength = Math.max(1000, n / 100);
    int[] cycle = new int[cycleMaxLength];
    int J = 0, i;
    S_k[0] = 2 * epsilon[0] - 1;
    for (i = 1; i < n; i++) {
        S_k[i] = S_k[i - 1] + 2 * epsilon[i] - 1;
        if (0 == S_k[i]) {
            J++;
            assertTrue("Exceeding the max number of cycles expected.", J <= cycleMaxLength);
            cycle[J] = i;
        }
    }
    if (0 != S_k[n - 1]) {
        J++;
        cycle[J] = n;
    }
    //int constraint = (int) Math.max(0.005 * Math.pow(n, 0.5), 500);
    int cycleStart = 0, cycleStop = cycle[1], j, k, b;
    double[][] nu = new double[6][8];
    int[] counter = { 0, 0, 0, 0, 0, 0, 0, 0 };
    double pi[][] = {
            { 0.0000000000, 0.00000000000, 0.00000000000, 0.00000000000, 0.00000000000, 0.0000000000 },
            { 0.5000000000, 0.25000000000, 0.12500000000, 0.06250000000, 0.03125000000, 0.0312500000 },
            { 0.7500000000, 0.06250000000, 0.04687500000, 0.03515625000, 0.02636718750, 0.0791015625 },
            { 0.8333333333, 0.02777777778, 0.02314814815, 0.01929012346, 0.01607510288, 0.0803755143 },
            { 0.8750000000, 0.01562500000, 0.01367187500, 0.01196289063, 0.01046752930, 0.0732727051 } };
    for (k = 0; k < 6; k++) {
        for (i = 0; i < 8; i++) {
            nu[k][i] = 0;
        }
    }
    for (j = 1; j <= J; j++) {
        for (i = 0; i < 8; i++) {
            counter[i] = 0;
        }
        for (i = cycleStart; i < cycleStop; i++) {
            if ((S_k[i] >= 1 && S_k[i] <= 4) || (S_k[i] >= -4 && S_k[i] <= -1)) {
                if (S_k[i] < 0)
                    b = 4;
                else
                    b = 3;
                counter[S_k[i] + b]++;
            }
        }
        cycleStart = cycle[j] + 1;
        if (j < J)
            cycleStop = cycle[j + 1];

        for (i = 0; i < 8; i++) {
            if ((counter[i] >= 0) && (counter[i] <= 4))
                nu[counter[i]][i]++;
            else if (counter[i] >= 5)
                nu[5][i]++;
        }
    }
    int x;
    double sum, p_value;
    for (i = 0; i < 8; i++) {
        x = stateX[i];
        sum = 0.;
        for (k = 0; k < 6; k++)
            sum += Math.pow(nu[k][i] - J * pi[(int) Math.abs(x)][k], 2) / (J * pi[(int) Math.abs(x)][k]);
        p_value = Gamma.regularizedGammaQ(2.5, sum / 2.0);
        assertTrue("RNG test failed, test random excursions.", p_value >= 0.01);
    }
}