using MathNet.Numerics.Distributions;

namespace IStation.Test
{
    public static class CorrelationHelper
    {

        /// <summary>
        /// 计算皮尔逊相关系数及显著性p值
        /// </summary>
        public static (double r, double pValue) PearsonTest(double[] x, double[] y)
        {
            // 计算皮尔逊相关系数
            double r = Pearson(x, y);

            // 计算显著性p值（双尾t检验）
            int n = x.Length;
            double t = r * Math.Sqrt((n - 2) / (1 - r * r));
            var tDist = new StudentT(location: 0, scale: 1, freedom: n - 2);
            double pValue = 2 * (1 - tDist.CumulativeDistribution(Math.Abs(t)));

            return (r, pValue);
        }

        /// <summary>
        /// 计算斯皮尔曼等级相关系数及显著性p值
        /// </summary>
        public static (double rho, double pValue) SpearmanTest(double[] x, double[] y)
        {
            // 计算斯皮尔曼相关系数
            double rho = Spearman(x, y);

            // 大样本正态近似法计算p值
            int n = x.Length;
            double z = rho * Math.Sqrt(n - 1);
            var normalDist = new Normal();
            double pValue = 2 * (1 - normalDist.CumulativeDistribution(Math.Abs(z)));

            return (rho, pValue);
        }


        /// <summary>
        /// 计算皮尔逊积差相关系数
        /// </summary>
        /// <param name="dataA">数据样本A</param>
        /// <param name="dataB">数据样本B</param>
        /// <returns>皮尔逊积差相关系数</returns>
        public static double Pearson(IEnumerable<double> dataA, IEnumerable<double> dataB)
        {
            int n = 0;
            double r = 0.0;
            double meanA = 0;
            double meanB = 0;
            double varA = 0;
            double varB = 0;

            using (var ieA = dataA.GetEnumerator())
            using (var ieB = dataB.GetEnumerator())
            {
                while (ieA.MoveNext())
                {
                    if (!ieB.MoveNext())
                    {
                        throw new ArgumentOutOfRangeException("dataB", "数据长度不一致");
                    }

                    double currentA = ieA.Current;
                    double currentB = ieB.Current;

                    double deltaA = currentA - meanA;
                    double scaleDeltaA = deltaA / ++n;

                    double deltaB = currentB - meanB;
                    double scaleDeltaB = deltaB / n;

                    meanA += scaleDeltaA;
                    meanB += scaleDeltaB;

                    varA += scaleDeltaA * deltaA * (n - 1);
                    varB += scaleDeltaB * deltaB * (n - 1);
                    r += (deltaA * deltaB * (n - 1)) / n;
                }

                if (ieB.MoveNext())
                {
                    throw new ArgumentOutOfRangeException("dataA", "数据长度不一致");
                }
            }

            return r / Math.Sqrt(varA * varB);
        }

        /// <summary>
        /// 计算斯皮尔曼等级相关系数
        /// </summary>
        public static double Spearman(IEnumerable<double> dataA, IEnumerable<double> dataB)
        {
            return Pearson(Rank(dataA), Rank(dataB));
        }

        private static IEnumerable<double> Rank(IEnumerable<double> sequence)
        {
            var sorted = new List<double>(sequence);
            sorted.Sort();

            return (IEnumerable<double>)sequence.Select(x => (double)sorted.IndexOf(x) + 1);
        }
    }
}