/*

* Copyright (c) NM LTD.

* https://nm.dev/

*

* THIS SOFTWARE IS LICENSED, NOT SOLD.

*

* YOU MAY USE THIS SOFTWARE ONLY AS DESCRIBED IN THE LICENSE.

* IF YOU ARE NOT AWARE OF AND/OR DO NOT AGREE TO THE TERMS OF THE LICENSE,

* DO NOT USE THIS SOFTWARE.

*

* THE SOFTWARE IS PROVIDED "AS IS", WITH NO WARRANTY WHATSOEVER,

* EITHER EXPRESS OR IMPLIED, INCLUDING, WITHOUT LIMITATION,

* ANY WARRANTIES OF ACCURACY, ACCESSIBILITY, COMPLETENESS,

* FITNESS FOR A PARTICULAR PURPOSE, MERCHANTABILITY, NON-INFRINGEMENT,

* TITLE AND USEFULNESS.

*

* IN NO EVENT AND UNDER NO LEGAL THEORY,

* WHETHER IN ACTION, CONTRACT, NEGLIGENCE, TORT, OR OTHERWISE,

* SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR

* ANY CLAIMS, DAMAGES OR OTHER LIABILITIES,

* ARISING AS A RESULT OF USING OR OTHER DEALINGS IN THE SOFTWARE.

*/

package dev.nm.nmj;

import java.util.List;

import java.util.ArrayList;

import java.util.Arrays;

import dev.nm.algebra.linear.matrix.doubles.Matrix;

import dev.nm.algebra.linear.matrix.doubles.matrixtype.dense.DenseMatrix;

import dev.nm.algebra.linear.matrix.doubles.matrixtype.dense.triangle.SymmetricMatrix;

import dev.nm.algebra.linear.matrix.doubles.operation.Inverse;

import dev.nm.algebra.linear.vector.doubles.Vector;

import dev.nm.algebra.linear.vector.doubles.dense.DenseVector;

import dev.nm.number.DoubleUtils;

import dev.nm.stat.covariance.*;

import dev.nm.stat.covariance.covarianceselection.CovarianceSelectionProblem;

import dev.nm.stat.covariance.covarianceselection.lasso.CovarianceSelectionGLASSOFAST;

import dev.nm.stat.covariance.covarianceselection.lasso.CovarianceSelectionLASSO;

import dev.nm.stat.covariance.nlshrink.LedoitWolf2016;

import dev.nm.stat.descriptive.correlation.CorrelationMatrix;

import dev.nm.stat.descriptive.correlation.SpearmanRankCorrelation;

import dev.nm.stat.descriptive.covariance.Covariance;

import dev.nm.stat.descriptive.covariance.SampleCovariance;

import dev.nm.stat.descriptive.moment.Kurtosis;

import dev.nm.stat.descriptive.moment.Mean;

import dev.nm.stat.descriptive.moment.Moments;

import dev.nm.stat.descriptive.moment.Skewness;

import dev.nm.stat.descriptive.moment.Variance;

import dev.nm.stat.descriptive.moment.weighted.WeightedMean;

import dev.nm.stat.descriptive.moment.weighted.WeightedVariance;

import dev.nm.stat.descriptive.rank.Max;

import dev.nm.stat.descriptive.rank.Min;

import dev.nm.stat.descriptive.rank.Quantile;

import dev.nm.stat.descriptive.rank.Rank;

import dev.nm.stat.distribution.discrete.ProbabilityMassQuantile;

import dev.nm.stat.distribution.univariate.WeibullDistribution;

import dev.nm.stat.distribution.discrete.ProbabilityMassFunction.Mass;

import dev.nm.stat.distribution.multivariate.DirichletDistribution;

import dev.nm.stat.distribution.multivariate.MultinomialDistribution;

import dev.nm.stat.distribution.multivariate.MultivariateNormalDistribution;

import dev.nm.stat.distribution.multivariate.MultivariateProbabilityDistribution;

import dev.nm.stat.distribution.multivariate.MultivariateTDistribution;

import dev.nm.stat.distribution.univariate.BetaDistribution;

import dev.nm.stat.distribution.univariate.BinomialDistribution;

import dev.nm.stat.distribution.univariate.ChiSquareDistribution;

import dev.nm.stat.distribution.univariate.EmpiricalDistribution;

import dev.nm.stat.distribution.univariate.ExponentialDistribution;

import dev.nm.stat.distribution.univariate.FDistribution;

import dev.nm.stat.distribution.univariate.GammaDistribution;

import dev.nm.stat.distribution.univariate.LogNormalDistribution;

import dev.nm.stat.distribution.univariate.NormalDistribution;

import dev.nm.stat.distribution.univariate.PoissonDistribution;

import dev.nm.stat.distribution.univariate.RayleighDistribution;

import dev.nm.stat.distribution.univariate.TDistribution;

import dev.nm.stat.factor.factoranalysis.FAEstimator;

import dev.nm.stat.factor.factoranalysis.FactorAnalysis;

import dev.nm.stat.factor.pca.PCA;

import dev.nm.stat.factor.pca.PCAbyEigen;

import dev.nm.stat.factor.pca.PCAbySVD;

import dev.nm.stat.hmm.ForwardBackwardProcedure;

import dev.nm.stat.hmm.HmmInnovation;

import dev.nm.stat.hmm.Viterbi;

import dev.nm.stat.hmm.discrete.BaumWelch;

import dev.nm.stat.hmm.discrete.DiscreteHMM;

import dev.nm.stat.hmm.mixture.MixtureHMM;

import dev.nm.stat.hmm.mixture.MixtureHMMEM;

import dev.nm.stat.hmm.mixture.distribution.NormalMixtureDistribution;

import dev.nm.stat.markovchain.SimpleMC;

import dev.nm.stat.random.rng.univariate.normal.StandardNormalRNG;

import dev.nm.stat.test.distribution.AndersonDarling;

import dev.nm.stat.test.distribution.CramerVonMises2Samples;

import dev.nm.stat.test.distribution.kolmogorov.KolmogorovSmirnov;

import dev.nm.stat.test.distribution.kolmogorov.KolmogorovSmirnov1Sample;

import dev.nm.stat.test.distribution.kolmogorov.KolmogorovSmirnov2Samples;

import dev.nm.stat.test.distribution.normality.DAgostino;

import dev.nm.stat.test.distribution.normality.JarqueBera;

import dev.nm.stat.test.distribution.normality.Lilliefors;

import dev.nm.stat.test.distribution.normality.ShapiroWilk;

import dev.nm.stat.test.distribution.pearson.ChiSquareIndependenceTest;

import dev.nm.stat.test.mean.OneWayANOVA;

import dev.nm.stat.test.mean.T;

import dev.nm.stat.test.rank.KruskalWallis;

import dev.nm.stat.test.rank.SiegelTukey;

import dev.nm.stat.test.rank.VanDerWaerden;

import dev.nm.stat.test.rank.wilcoxon.WilcoxonSignedRank;

import static java.lang.Math.abs;

import static java.lang.Math.sqrt;

/**

* Numerical Methods Using Java: For Data Science, Analysis, and Engineering

*

* @author haksunli

* @see

* https://www.amazon.com/Numerical-Methods-Using-Java-Engineering/dp/1484267966

* https://nm.dev/

*/

public class Chapter12 {

public static void main(String[] args) throws Exception {

System.out.println("Chapter 12 demos");

Chapter12 chapter12 = new Chapter12();

chapter12.sample_statistics();

chapter12.rank();

chapter12.quantile();

chapter12.LedoitWolf2004();

chapter12.LedoitWolf2016();

chapter12.normal_distributions();

chapter12.lognormal_distributions();

chapter12.exponential_distribution();

chapter12.Poisson_distribution();

chapter12.binomial_distribution();

chapter12.t_distribution();

chapter12.F_distribution();

chapter12.chi_square_distribution();

chapter12.Rayleigh_distribution();

chapter12.gamma_distribution();

chapter12.beta_distribution();

chapter12.Weibull_distribution();

chapter12.empirical_distribution();

chapter12.multivariate_normal_distribution();

chapter12.multivariate_t_distribution();

chapter12.Dirichlet_distribution();

chapter12.multinomial_distribution();

chapter12.hypothesis_testing();

chapter12.Shapiro_Wilk_test();

chapter12.Jarque_Bera_test();

chapter12.DAgostino_test();

chapter12.Lilliefors_test();

chapter12.Kolmogorov_Smirnov_test();

chapter12.Anderson_Darling_test();

chapter12.Cramer_Von_Mises_test();

chapter12.Chi_square_independence_test();

chapter12.t_test();

chapter12.one_way_ANOVA();

chapter12.Kruskal_Wallis_test();

chapter12.Wilcoxon_signed_rank_test();

chapter12.Siegel_Tukey_test();

chapter12.Van_Der_Waerden_test();

chapter12.DTMC();

chapter12.HMM1();

chapter12.HMM2();

chapter12.HMM3();

chapter12.PCA_engen();

chapter12.PCA_svd();

chapter12.factor_analysis();

chapter12.covariance_selection_LASSO();

}

public void covariance_selection_LASSO() {

System.out.println("covariance selection using LASSO");

// generate random samples from standard normal distribution

StandardNormalRNG rnorm = new StandardNormalRNG();

rnorm.seed(1234567890L);

int nRows = 50;

int nCols = 10;

Matrix X = new DenseMatrix(nRows, nCols);

for (int i = 1; i <= nRows; ++i) {

for (int j = 1; j <= nCols; ++j) {

X.set(i, j, rnorm.nextDouble());

}

}

// sample covariance matrix

Matrix S = new SampleCovariance(X);

System.out.println("sample covariance:");

System.out.println(S);

Matrix S_inv = new Inverse(S);

System.out.println("inverse sample covariance:");

System.out.println(S_inv);

// the penalty parameter

double rho = 0.03;

CovarianceSelectionProblem problem

= new CovarianceSelectionProblem(S, rho);

long time1 = System.currentTimeMillis();

CovarianceSelectionLASSO lasso

= new CovarianceSelectionLASSO(problem, 1e-5);

Matrix sigma = lasso.covariance();

time1 = System.currentTimeMillis() - time1;

System.out.println("estimated sigma:");

System.out.println(sigma);

System.out.println("inverse sigma:");

Matrix sigma_inv = lasso.inverseCovariance();

System.out.println(sigma_inv);

long time2 = System.currentTimeMillis();

CovarianceSelectionGLASSOFAST lasso2 = new CovarianceSelectionGLASSOFAST(problem);

Matrix sigma2 = lasso2.covariance();

time2 = System.currentTimeMillis() - time2;

System.out.println("CovarianceSelectionLASSO took " + time1 + " millisecs");

System.out.println("CovarianceSelectionGLASSOFAST took " + time2 + " millisecs");

}

// data set from R

private static final Matrix R_data

= new DenseMatrix(new double[][]{

{1., 1., 3., 3., 1., 1.},

{1., 2., 3., 3., 1., 1.},

{1., 1., 3., 4., 1., 1.},

{1., 1., 3., 3., 1., 2.},

{1., 1., 3., 3., 1., 1.},

{1., 1., 1., 1., 3., 3.},

{1., 2., 1., 1., 3., 3.},

{1., 1., 1., 2., 3., 3.},

{1., 2., 1., 1., 3., 4.},

{1., 1., 1., 1., 3., 3.},

{3., 3., 1., 1., 1., 1.},

{3., 4., 1., 1., 1., 1.},

{3., 3., 1., 2., 1., 1.},

{3., 3., 1., 1., 1., 2.},

{3., 3., 1., 1., 1., 1.},

{4., 4., 5., 5., 6., 6.},

{5., 6., 4., 6., 4., 5.},

{6., 5., 6., 4., 5., 4.}

});

public void factor_analysis() {

System.out.println("factor analysis");

// number of hidden factors

int nFactors = 3;

FactorAnalysis factor_analysis

= new FactorAnalysis(

R_data,

nFactors,

FactorAnalysis.ScoringRule.THOMSON // specify the scoring rule

);

System.out.println("number of observations = " + factor_analysis.nObs());

System.out.println("number of variables = " + factor_analysis.nVariables());

System.out.println("number of factors = " + factor_analysis.nFactors());

// covariance matrix

Matrix S = factor_analysis.S();

System.out.println("covariance matrix:");

System.out.println(S);

FAEstimator estimators = factor_analysis.getEstimators(700);

double fitted = estimators.logLikelihood();

System.out.println("log-likelihood of the fitting = " + fitted);

Vector uniqueness = estimators.psi();

System.out.println("uniqueness = " + uniqueness);

int dof = estimators.dof();

System.out.println("degree of freedom = " + dof);

// the factor loadings

Matrix loadings = estimators.loadings();

System.out.println("factor loadings:");

System.out.println(loadings);

double testStats = estimators.statistics();

System.out.println("test statistics = " + testStats);

double pValue = estimators.pValue();

System.out.println("p-value = " + pValue);

Matrix scores = estimators.scores();

System.out.println("scores:");

System.out.println(scores);

}

// this is the US arrest data from R

private static final DenseMatrix USArrests

= new DenseMatrix(new double[][]{

{13.2, 236, 58, 21.2},

{10.0, 263, 48, 44.5},

{8.1, 294, 80, 31.0},

{8.8, 190, 50, 19.5},

{9.0, 276, 91, 40.6},

{7.9, 204, 78, 38.7},

{3.3, 110, 77, 11.1},

{5.9, 238, 72, 15.8},

{15.4, 335, 80, 31.9},

{17.4, 211, 60, 25.8},

{5.3, 46, 83, 20.2},

{2.6, 120, 54, 14.2},

{10.4, 249, 83, 24.0},

{7.2, 113, 65, 21.0},

{2.2, 56, 57, 11.3},

{6.0, 115, 66, 18.0},

{9.7, 109, 52, 16.3},

{15.4, 249, 66, 22.2},

{2.1, 83, 51, 7.8},

{11.3, 300, 67, 27.8},

{4.4, 149, 85, 16.3},

{12.1, 255, 74, 35.1},

{2.7, 72, 66, 14.9},

{16.1, 259, 44, 17.1},

{9.0, 178, 70, 28.2},

{6.0, 109, 53, 16.4},

{4.3, 102, 62, 16.5},

{12.2, 252, 81, 46.0},

{2.1, 57, 56, 9.5},

{7.4, 159, 89, 18.8},

{11.4, 285, 70, 32.1},

{11.1, 254, 86, 26.1},

{13.0, 337, 45, 16.1},

{0.8, 45, 44, 7.3},

{7.3, 120, 75, 21.4},

{6.6, 151, 68, 20.0},

{4.9, 159, 67, 29.3},

{6.3, 106, 72, 14.9},

{3.4, 174, 87, 8.3},

{14.4, 279, 48, 22.5},

{3.8, 86, 45, 12.8},

{13.2, 188, 59, 26.9},

{12.7, 201, 80, 25.5},

{3.2, 120, 80, 22.9},

{2.2, 48, 32, 11.2},

{8.5, 156, 63, 20.7},

{4.0, 145, 73, 26.2},

{5.7, 81, 39, 9.3},

{2.6, 53, 66, 10.8},

{6.8, 161, 60, 15.6}

});

public void PCA_svd() {

System.out.println("PCA by SVD");

// run PCA on the data using SVD

PCA pca = new PCAbySVD(

USArrests // the data set in matrix form

);

// number of factors

int p = pca.nFactors();

// number of observations

int n = pca.nObs();

Vector mean = pca.mean();

Vector scale = pca.scale();

Vector sdev = pca.sdPrincipalComponents();

Matrix loadings = pca.loadings();

Vector proportion = pca.proportionVar();

Vector cumprop = pca.cumulativeProportionVar();

System.out.println("number of factors = " + p);

System.out.println("number of observations = " + n);

System.out.println("mean: " + mean);

System.out.println("scale: " + scale);

// The standard deviations differ by a factor of sqrt(50 / 49),

// since we use divisor (nObs - 1) for the sample covariance matrix

System.out.println("standard deviation: " + sdev);

// The signs of the columns of the loading are arbitrary.

System.out.println("loading: ");

System.out.println(loadings);

// the proportion of variance in each dimension

System.out.println("proportion of variance: " + proportion);

System.out.println("cumulative proportion of variance: " + cumprop);

}

public void PCA_engen() {

System.out.println("PCA by eigen decomposition");

// run PCA on the data using eigen decomposition

PCA pca = new PCAbyEigen(

USArrests, // the data set in matrix form

false // use covariance matrix instead of correlation matrix

);

// number of factors

int p = pca.nFactors();

// number of observations

int n = pca.nObs();

Vector mean = pca.mean();

Vector scale = pca.scale();

Vector sdev = pca.sdPrincipalComponents();

Matrix loadings = pca.loadings();

Vector proportion = pca.proportionVar();

Vector cumprop = pca.cumulativeProportionVar();

Matrix scores = pca.scores();

System.out.println("number of factors = " + p);

System.out.println("number of observations = " + n);

System.out.println("mean: " + mean);

System.out.println("scale: " + scale);

// The standard deviations differ by a factor of sqrt(50 / 49),

// since we use divisor (nObs - 1) for the sample covariance matrix

System.out.println("standard deviation: " + sdev);

// The signs of the columns of the loading are arbitrary.

System.out.println("loading: ");

System.out.println(loadings);

// the proportion of variance in each dimension

System.out.println("proportion of variance: " + proportion);

System.out.println("cumulative proportion of variance: " + cumprop);

// System.out.println("score: ");

// System.out.println(scores);

}

public void HMM3() {

System.out.println("learning hidden Markov model with normal distribution");

// the initial probabilities

Vector PI0 = new DenseVector(new double[]{0., 0., 1.});

// the transition probabilities

Matrix A0 = new DenseMatrix(new double[][]{

{0.4, 0.2, 0.4},

{0.3, 0.2, 0.5},

{0.25, 0.25, 0.5}

});

// the conditional normal distributions

NormalMixtureDistribution.Lambda[] lambda0 = new NormalMixtureDistribution.Lambda[]{

new NormalMixtureDistribution.Lambda(0., .5), // (mu, sigma)

new NormalMixtureDistribution.Lambda(0., 1.), // medium volatility

new NormalMixtureDistribution.Lambda(0., 2.5) // high volatility

};

// the original HMM: a model of daily stock returns in 3 regimes

MixtureHMM model0 = new MixtureHMM(PI0, A0, new NormalMixtureDistribution(lambda0));

model0.seed(1234567890L);

// generate a sequence of observations from the HMM

int T = 10000;

HmmInnovation[] innovations = new HmmInnovation[T];

double[] observations = new double[T];

for (int t = 0; t < T; ++t) {

innovations[t] = model0.next();

observations[t] = innovations[t].observation();

}

System.out.println("observations: ");

for (int t = 1; t <= 100; ++t) {

System.out.print(observations[t] + ", ");

if (t % 20 == 0) {

System.out.println("");

}

}

System.out.println("");

// learn an HMM from the observations

MixtureHMM model1

= new MixtureHMMEM(observations, model0, 1e-5, 20); // using true parameters as initial estimates

Matrix A1 = model1.A();

NormalMixtureDistribution.Lambda[] lambda1

= ((NormalMixtureDistribution) model1.getDistribution()).getParams();

System.out.println("original transition probabilities");

System.out.println(A0);

System.out.println("learned transition probabilities");

System.out.println(A1);

for (int i = 0; i < lambda0.length; ++i) {

System.out.println(String.format("compare mu: %f vs %f", lambda0[i].mu, lambda1[i].mu));

System.out.println(String.format("compare sigma: %f vs %f", lambda0[i].sigma, lambda1[i].sigma));

}

}

public void HMM2() {

System.out.println("learning hidden Markov model");

// generate a sequence of observations from a HMM

// the initial probabilities for 2 states

DenseVector PI = new DenseVector(

new double[]{0.6, 0.4}

);

// the transition probabilities

DenseMatrix A = new DenseMatrix(new double[][]{

{0.7, 0.3},

{0.4, 0.6}

});

// the observation probabilities; 3 possible outcomes for 2 states

DenseMatrix B = new DenseMatrix(new double[][]{

{0.5, 0.4, 0.1},

{0.1, 0.3, 0.6}

});

// construct an HMM1

DiscreteHMM model = new DiscreteHMM(PI, A, B);

model.seed(1234507890L, 1234507891L);

// generate the observations

int T = 10000;

HmmInnovation[] innovations = new HmmInnovation[T];

int[] states = new int[T];

int[] observations = new int[T];

for (int t = 0; t < T; ++t) {

innovations[t] = model.next();

states[t] = innovations[t].state();

observations[t] = (int) innovations[t].observation();

}

System.out.println("observations: ");

for (int t = 1; t <= 100; ++t) {

System.out.print(observations[t] + ", ");

if (t % 20 == 0) {

System.out.println("");

}

}

System.out.println("");

// learn the HMM from observations

DenseVector PI_0 = new DenseVector(new double[]{0.5, 0.5}); // initial guesses

DenseMatrix A_0 = new DenseMatrix(new double[][]{ // initial guesses

{0.5, 0.5},

{0.5, 0.5}

});

DenseMatrix B_0 = new DenseMatrix(new double[][]{ // initial guesses

{0.60, 0.20, 0.20},

{0.20, 0.20, 0.60}

});

DiscreteHMM model_0 = new DiscreteHMM(PI_0, A_0, B_0); // initial guesses

// training

int nIterations = 40;

for (int i = 1; i <= nIterations; ++i) {

model_0 = BaumWelch.train(observations, model_0);

}

// training results

System.out.println("estimated transition probabilities: ");

System.out.println(model_0.A()); // should be close to A

System.out.println("(observation) conditional probabilities: ");

System.out.println(model_0.B()); // should be close to B

}

public void HMM1() {

System.out.println("hidden Markov model");

// the initial probabilities for 2 states

DenseVector PI = new DenseVector(

new double[]{0.6, 0.4}

);

// the transition probabilities

DenseMatrix A = new DenseMatrix(new double[][]{

{0.7, 0.3},

{0.4, 0.6}

});

// the observation probabilities; 3 possible outcomes for 2 states

DenseMatrix B = new DenseMatrix(new double[][]{

{0.5, 0.4, 0.1},

{0.1, 0.3, 0.6}

});

// construct an HMM1

DiscreteHMM hmm = new DiscreteHMM(PI, A, B);

// the realized observations

double[] observations = new double[]{1, 2, 3};

// run the forward-backward algorithm

ForwardBackwardProcedure fb = new ForwardBackwardProcedure(hmm, observations);

for (int t = 1; t <= observations.length; ++t) {

System.out.println(String.format(

"the *scaled* forward probability, alpha, in each state at time %d: %s",

t,

fb.scaledAlpha(t)

));

}

// run the Viterbi algorithm to find the most likely sequence of hidden states

Viterbi viterbi = new Viterbi(hmm);

int[] viterbi_states = viterbi.getViterbiStates(observations);

System.out.println("the Viterbi states: " + Arrays.toString(viterbi_states));

}

public void DTMC() {

System.out.println("discrete time Markov chain");

// the stochastic matrix of transition probabilities

Matrix A = new DenseMatrix(new double[][]{

{0.4, 0.2, 0.4},

{0.3, 0.2, 0.5},

{0.25, 0.25, 0.5}

});

// start in state 3

Vector I = new DenseVector(0., 0., 1.);

SimpleMC MC = new SimpleMC(I, A);

Vector PI = SimpleMC.getStationaryProbabilities(A);

System.out.println("the stationary distribution = " + PI);

// simulate the next 9 steps

System.out.println("time 0 = " + 3);

for (int i = 1; i < 10; ++i) {

int state = MC.nextState();

System.out.println(String.format("time %d = %d", i, state));

}

}

public void Van_Der_Waerden_test() {

System.out.println("Van der Waerden test");

double[][] samples = new double[3][];

samples[0] = new double[]{8, 10, 9, 10, 9};

samples[1] = new double[]{7, 8, 5, 8, 5};

samples[2] = new double[]{4, 8, 7, 5, 7};

VanDerWaerden test = new VanDerWaerden(samples);

System.out.println("H0: " + test.getNullHypothesis());

System.out.println("H1: " + test.getAlternativeHypothesis());

System.out.println("test statistics = " + test.statistics());

System.out.println("p-value = " + test.pValue());

System.out.println("is null rejected at 5% = " + test.isNullRejected(0.05));

}

public void Siegel_Tukey_test() {

System.out.println("Siegel Tukey test");

double[] sample1 = new double[]{4, 16, 48, 51, 66, 98};

double[] sample2 = new double[]{33, 62, 84, 85, 88, 93, 97};

SiegelTukey test = new SiegelTukey(

sample1,

sample2,

0, // the hypothetical mean difference

true // use the exact Wilcoxon Rank Sum distribution rather than normal distribution

);

System.out.println("H0: " + test.getNullHypothesis());

System.out.println("H1: " + test.getAlternativeHypothesis());

System.out.println("test statistics = " + test.statistics());

System.out.println("p-value = " + test.pValue());

System.out.println("p-value, right sided = " + test.rightOneSidedPvalue());

System.out.println("p-value, left sided = " + test.leftOneSidedPvalue());

System.out.println("is null rejected at 5% = " + test.isNullRejected(0.05));

}

public void Wilcoxon_signed_rank_test() {

System.out.println("Wilcoxon signed rank test");

double[] sample1 = new double[]{1.3, 5.4, 7.6, 7.2, 3.5};

double[] sample2 = new double[]{2.7, 5.2, 6.3, 4.4, 9.8};

WilcoxonSignedRank test = new WilcoxonSignedRank(

sample1, sample2,

2, // the hypothetical median that the distribution is symmetric about

true // use the exact Wilcoxon rank sum distribution rather than normal distribution

);

System.out.println("H0: " + test.getNullHypothesis());

System.out.println("H1: " + test.getAlternativeHypothesis());

System.out.println("test statistics = " + test.statistics());

System.out.println("p-value = " + test.pValue());

System.out.println("p-value, right sided = " + test.rightOneSidedPvalue());

System.out.println("p-value, left sided = " + test.leftOneSidedPvalue());

System.out.println("is null rejected at 5% = " + test.isNullRejected(0.05));

}

public void Kruskal_Wallis_test() {

System.out.println("Kruskal Wallis test");

double[][] samples = new double[4][];

samples[0] = new double[]{1, 1, 7.6, 7.2, 3.5};

samples[1] = new double[]{2, 2, 6.3, 4.4, 9.8, 10.24};

samples[2] = new double[]{-9, -9, -4.33, -5.4};

samples[3] = new double[]{0.21, 0.21, 0.21, 0.86, 0.902, 0.663};

KruskalWallis test = new KruskalWallis(samples);

System.out.println("H0: " + test.getNullHypothesis());

System.out.println("H1: " + test.getAlternativeHypothesis());

System.out.println("test statistics = " + test.statistics());

System.out.println("p-value = " + test.pValue());

System.out.println("is null rejected at 5% = " + test.isNullRejected(0.05));

}

public void one_way_ANOVA() {

System.out.println("One-way ANOVA");

double[][] samples = new double[4][];

samples[0] = new double[]{1.3, 5.4, 7.6, 7.2, 3.5};

samples[1] = new double[]{2.7, 5.21, 6.3, 4.4, 9.8, 10.24};

samples[2] = new double[]{-2.3, -5.3, -4.33, -5.4};

samples[3] = new double[]{0.21, 0.34, 0.27, 0.86, 0.902, 0.663};

OneWayANOVA test = new OneWayANOVA(samples);

System.out.println("H0: " + test.getNullHypothesis());

System.out.println("H1: " + test.getAlternativeHypothesis());

System.out.println("test statistics = " + test.statistics());

System.out.println("p-value = " + test.pValue());

System.out.println("is null rejected at 5% = " + test.isNullRejected(0.05));

}

public void t_test() {

System.out.println("t test");

// the t-test

T test1 = new T(

new double[]{1, 3, 5, 2, 3, 5},

new double[]{2, 5, 6, 4, 9, 8},

true, // assume variances are equal

4 // the hypothetical mean-difference = 4 in the null hypothesis

);

System.out.println("H0: " + test1.getNullHypothesis());

System.out.println("H1: " + test1.getAlternativeHypothesis());

System.out.println("test statistics = " + test1.statistics());

System.out.println("1st mean = " + test1.mean1());

System.out.println("2nd mean = " + test1.mean2());

System.out.println("p-value = " + test1.pValue());

System.out.println("p-value, right sided = " + test1.rightOneSidedPvalue());

System.out.println("p-value, left sided = " + test1.leftOneSidedPvalue());

System.out.println(String.format("95%% confidence interval = (%f, %f)", test1.leftConfidenceInterval(0.95), test1.rightConfidenceInterval(0.95)));

System.out.println("97.5%% confidence interval = " + Arrays.toString(test1.confidenceInterval(0.975)));

System.out.println("is null rejected at 5% = " + test1.isNullRejected(0.05));

// Welch's t-test

T test2 = new T(

new double[]{1, 3, 5, 2, 3, 5},

new double[]{2, 5, 6, 4, 9, 8},

false, // assume variances are different

4 // the hypothetical mean-difference = 4 in the null hypothesis

);

System.out.println("test statistics = " + test2.statistics());

System.out.println("p-value = " + test2.pValue());

System.out.println("p-value, right sided = " + test2.rightOneSidedPvalue());

System.out.println("p-value, left sided = " + test2.leftOneSidedPvalue());

System.out.println(String.format("95%% confidence interval = (%f, %f)", test2.leftConfidenceInterval(0.95), test2.rightConfidenceInterval(0.95)));

System.out.println("97.5%% confidence interval = " + Arrays.toString(test2.confidenceInterval(0.975)));

System.out.println("is null rejected at 5% = " + test2.isNullRejected(0.05));

}

public void Chi_square_independence_test() {

System.out.println("Chi-square independence test");

// the attendance/absence vs. pass/fail counts

Matrix counts = new DenseMatrix(new double[][]{

{25, 6},

{8, 15}

});

ChiSquareIndependenceTest test1

= new ChiSquareIndependenceTest(

counts,

0,

// the asymptotic distribution is the Chi-square distribution

ChiSquareIndependenceTest.Type.ASYMPTOTIC

);

Matrix expected = ChiSquareIndependenceTest.getExpectedContingencyTable(

new int[]{31, 23}, // row sums

new int[]{33, 21} // column sums

);

System.out.println("the expected frequencies:");

System.out.println(expected);

System.out.println("H0: " + test1.getNullHypothesis());

System.out.println("H1: " + test1.getAlternativeHypothesis());

System.out.println("test statistics = " + test1.statistics());

System.out.println("p-value = " + test1.pValue());

System.out.println("is null rejected at 5% = " + test1.isNullRejected(0.05));

ChiSquareIndependenceTest test2

= new ChiSquareIndependenceTest(

counts,

100000,// number of simulation to compute the Fisher exact distribution

ChiSquareIndependenceTest.Type.EXACT // use the Fisher exact distribution

);

System.out.println("p-value = " + test2.pValue());

System.out.println("is null rejected at 5% = " + test2.isNullRejected(0.05));

}

public void Cramer_Von_Mises_test() {

System.out.println("Cramer Von Mises test");

// the samples

double[] x1 = new double[]{-0.54289848, 0.08999578, -1.77719573, -0.67991860, -0.65741590, -0.25776164, 1.02024626, 1.26434300, 0.51068476, -0.23998229};

double[] x2 = new double[]{1.7053818, 1.0260726, 1.7695157, 1.5650577, 1.4945107, 1.8593791, 2.1760302, -0.9728721, 1.4208313, 1.5892663};

CramerVonMises2Samples test = new CramerVonMises2Samples(x1, x2);

System.out.println("H0: " + test.getNullHypothesis());

System.out.println("H1: " + test.getAlternativeHypothesis());

System.out.println("test statistics = " + test.statistics());

System.out.println("p-value = " + test.pValue());

System.out.println("is null rejected at 5% = " + test.isNullRejected(0.05));

}

public void Anderson_Darling_test() {

System.out.println("Anderson Darling test");

// the samples

double[] x1 = new double[]{38.7, 41.5, 43.8, 44.5, 45.5, 46.0, 47.7, 58.0};

double[] x2 = new double[]{39.2, 39.3, 39.7, 41.4, 41.8, 42.9, 43.3, 45.8};

double[] x3 = new double[]{34.0, 35.0, 39.0, 40.0, 43.0, 43.0, 44.0, 45.0};

double[] x4 = new double[]{34.0, 34.8, 34.8, 35.4, 37.2, 37.8, 41.2, 42.8};

AndersonDarling test = new AndersonDarling(x1, x2, x3, x4);

System.out.println("H0: " + test.getNullHypothesis());

System.out.println("H1: " + test.getAlternativeHypothesis());

System.out.println("test statistics = " + test.statistics());

System.out.println("p-value = " + test.pValue());

System.out.println("alternative test statistics = " + test.statisticsAlternative());

System.out.println("alternative p-value = " + test.pValueAlternative());

System.out.println("is null rejected at 5% = " + test.isNullRejected(0.05));

}

public void Kolmogorov_Smirnov_test() {

System.out.println("Kolmogorov Smirnov test");

// one-sample KS test

KolmogorovSmirnov1Sample test1 = new KolmogorovSmirnov1Sample(

new double[]{ // with duplicates

1.2142038235675114, 0.8271665834857130, -2.2786245743283295, 0.8414895245471727,

-1.4327682855296735, -0.2501807766164897, -1.9512765152306415, 0.6963626117638846,

0.4741320101265005, 1.2142038235675114

},

new NormalDistribution(),

KolmogorovSmirnov.Side.TWO_SIDED // options are: TWO_SIDED, GREATER, LESS

);

System.out.println("H0: " + test1.getNullHypothesis());

System.out.println("test statistics = " + test1.statistics());

System.out.println("p-value = " + test1.pValue());

System.out.println("is null rejected at 5% = " + test1.isNullRejected(0.05));

// two-sample KS test

KolmogorovSmirnov2Samples test2 = new KolmogorovSmirnov2Samples(

new double[]{ // x = rnorm(10)

1.2142038235675114, 0.8271665834857130, -2.2786245743283295, 0.8414895245471727,

-1.4327682855296735, -0.2501807766164897, -1.9512765152306415, 0.6963626117638846,

0.4741320101265005, -1.2340784297133520

},

new double[]{ // x = rnorm(15)

1.7996197748754565, -1.1371109188816089, 0.8179707525071304, 0.3809791236763478,

0.1644848304811257, 0.3397412780581336, -2.2571685407244795, 0.4137315314876659,

0.7318687611171864, 0.9905218801425318, -0.4748590846019594, 0.8882674167954235,

1.0534065683777052, 0.2553123235884622, -2.3172807717538038},

KolmogorovSmirnov.Side.GREATER // options are: TWO_SIDED, GREATER, LESS

);

System.out.println("H0: " + test2.getNullHypothesis());

System.out.println("test statistics = " + test2.statistics());

System.out.println("p-value = " + test2.pValue());

System.out.println("is null rejected at 5% = " + test2.isNullRejected(0.05));

}

public void Lilliefors_test() {

System.out.println("Lilliefors test");

double[] sample = new double[]{-1.7, -1, -1, -.73, -.61, -.5, -.24, .45, .62, .81, 1, 5};

Lilliefors test = new Lilliefors(sample);

System.out.println("H0: " + test.getNullHypothesis());

System.out.println("H1: " + test.getAlternativeHypothesis());

System.out.println("test statistics = " + test.statistics());

System.out.println("p-value = " + test.pValue());

System.out.println("is null rejected at 5% = " + test.isNullRejected(0.05));

}

public void DAgostino_test() {

System.out.println("D'Agostino's test");

double[] samples = new double[]{

39, 35, 33, 33, 32, 30, 30, 30, 28, 28,

27, 27, 27, 27, 27, 26, 26, 26, 26, 26,

26, 25, 25, 25, 25, 25, 25, 24, 24, 24,

24, 24, 23, 23, 23, 23, 23, 23, 23, 23,

23, 23, 23, 23, 23, 22, 22, 22, 22, 21,

21, 21, 21, 21, 21, 21, 20, 20, 19, 19,

18, 16

};

DAgostino test = new DAgostino(samples);

System.out.println("H0: " + test.getNullHypothesis());

System.out.println("H1: " + test.getAlternativeHypothesis());

System.out.println("skewness test statistics " + test.Z1());

System.out.println("p-value for skewness test = " + test.pvalueZ1());

System.out.println("kurtosis test statistics " + test.Z2());

System.out.println("test statistics = " + test.statistics());

System.out.println("p-value = " + test.pValue());

System.out.println("is null rejected at 5% = " + test.isNullRejected(0.05));

}

public void Jarque_Bera_test() {

System.out.println("Jarque-Bera test");

double[] samples = new double[]{

39, 35, 33, 33, 32, 30, 30, 30, 28, 28,

27, 27, 27, 27, 27, 26, 26, 26, 26, 26,

26, 25, 25, 25, 25, 25, 25, 24, 24, 24,

24, 24, 23, 23, 23, 23, 23, 23, 23, 23,

23, 23, 23, 23, 23, 22, 22, 22, 22, 21,

21, 21, 21, 21, 21, 21, 20, 20, 19, 19,

18, 16

};

JarqueBera test = new JarqueBera(

samples,

false // not using the exact Jarque-Bera distribution

);

System.out.println("H0: " + test.getNullHypothesis());

System.out.println("H1: " + test.getAlternativeHypothesis());

System.out.println("test statistics = " + test.statistics());

System.out.println("p-value = " + test.pValue());

System.out.println("is null rejected at 5% = " + test.isNullRejected(0.05));

}

public void Shapiro_Wilk_test() {

System.out.println("Shapiro-Wilk test");

double[] sample = new double[]{-1.7, -1, -1, -.73, -.61, -.5, -.24, .45, .62, .81, 1, 5};

ShapiroWilk test = new ShapiroWilk(sample);

System.out.println("H0: " + test.getNullHypothesis());

System.out.println("H1: " + test.getAlternativeHypothesis());

System.out.println("test statistics = " + test.statistics());

System.out.println("p-value = " + test.pValue());

System.out.println("is null rejected at 5% = " + test.isNullRejected(0.05));

}

public void hypothesis_testing() {

System.out.println("hypothesis testing");

int n = 100;

BinomialDistribution dist2 = new BinomialDistribution(

n,

0.5 // p

);

double stdev = sqrt(dist2.variance()) / n;

System.out.println("standard deviation = " + stdev);

double z_score = (0.37 - 0.5) / stdev;

System.out.println("z-score = " + z_score);

double p_value = new NormalDistribution() // default ctor for standard normal distribution

.cdf(z_score);

System.out.println("p-value = " + p_value);

}

public void multinomial_distribution() {

System.out.println("multinomial distribution");

// k = 3, each of the 3 probabilities of success

double[] prob = new double[]{0.1, 0.2, 0.7};

int n = 100;

MultinomialDistribution dist

= new MultinomialDistribution(n, prob);

// an outcome of the n trials

Vector x = new DenseVector(new double[]{10, 20, 70});

System.out.println(String.format("f(%s) = %f", x, dist.density(x)));

}

public void Dirichlet_distribution() {

System.out.println("Dirichlet distribution");

// the parameters

double[] a = new double[]{1, 2, 3, 4, 5};

DirichletDistribution dist = new DirichletDistribution(a);

Vector x = new DenseVector(0.1, 0.2, 0.3, 0.2, 0.2);

System.out.println(String.format("f(%s) = %f", x, dist.density(x)));

}

public void multivariate_t_distribution() {

System.out.println("multivariate t distribution");

int p = 2; // dimension

Vector mu = new DenseVector(1., 2.); // mean

Matrix Sigma = new DenseMatrix(p, p).ONE(); // scale matrix

int v = 1; // degree of freedom

MultivariateTDistribution t

= new MultivariateTDistribution(v, mu, Sigma);

Vector x = new DenseVector(1.23, 4.56);

System.out.println(String.format("f(%s) = %f", x, t.density(x)));

v = 2;

t = new MultivariateTDistribution(v, mu, Sigma);

x = new DenseVector(1.23, 4.56);

System.out.println(String.format("f(%s) = %f", x, t.density(x)));

v = 3;

t = new MultivariateTDistribution(v, mu, Sigma);

x = new DenseVector(1.23, 4.56);

System.out.println(String.format("f(%s) = %f", x, t.density(x)));

v = 4;

t = new MultivariateTDistribution(v, mu, Sigma);

x = new DenseVector(1.23, 4.56);

System.out.println(String.format("f(%s) = %f", x, t.density(x)));

v = 5;

t = new MultivariateTDistribution(v, mu, Sigma);

x = new DenseVector(1.23, 4.56);

System.out.println(String.format("f(%s) = %f", x, t.density(x)));

v = 6;

t = new MultivariateTDistribution(v, mu, Sigma);

x = new DenseVector(1.23, 4.56);

System.out.println(String.format("f(%s) = %f", x, t.density(x)));

}

public void multivariate_normal_distribution() {