C#，数值计算——用于从连续的数据值流估计任意分位数的计算方法与源程序

1 分位数Quantile

分位数（Quantile），亦称分位点，是指将一个随机变量的概率分布范围分为几个等份的数值点，常用的有中位数（即二分位数）、四分位数、百分位数等。

2 常见各类分位数

2.1 二分位数

对于有限的数集，可以通过把所有观察值高低排序后找出正中间的一个作为中位数。如果观察值有偶数个，则中位数不唯一，通常取最中间的两个数值的平均数作为中位数，即二分位数。
一个数集中最多有一半的数值小于中位数，也最多有一半的数值大于中位数。如果大于和小于中位数的数值个数均少于一半，那么数集中必有若干值等同于中位数。
计算有限个数的数据的二分位数的方法是：把所有的同类数据按照大小的顺序排列。如果数据的个数是奇数，则中间那个数据就是这群数据的中位数；如果数据的个数是偶数，则中间那2个数据的算术平均值就是这群数据的中位数。

2.2 四分位数

四分位数（Quartile）是统计学中分位数的一种，即把所有数值由小到大排列并分成四等份，处于三个分割点位置的数值就是四分位数。
1）第一四分位数(Q1)，又称“较小四分位数”，等于该样本中所有数值由小到大排列后第25%的数字；
2）第二四分位数(Q2)，又称“中位数”，等于该样本中所有数值由小到大排列后第50%的数字；
3）第三四分位数(Q3)，又称“较大四分位数”，等于该样本中所有数值由小到大排列后第75%的数字。
第三四分位数与第一四分位数的差距又称四分位距。

2.3 百分位数

百分位数，统计学术语，如果将一组数据从小到大排序，并计算相应的累计百分位，则某一百分位所对应数据的值就称为这一百分位的百分位数。运用在教育统计学中，例如表现测验成绩时，称PR值。

3 分位数的应用

分位数回归思想的提出至今已经有近30多年了，经过这近30多年的发展，分位数回归在理论和方法上都越来越成熟，并被广泛应用于多种学科中。它对于实际问题能提供更加全面的分析,无论是线性模型还是非线性模型，分位数回归都是一种很好的工具，它对一般回归模型做了有益的补充。

分位数回归是对以古典条件均值模型为基础的最小二乘法的延伸，它用几个分位函数来估计整体模型。分位数回归法的特殊情况就是中位数回归(最小一乘回归)，用对称权重解决残差最小化问题，而其他条件分位数回归则需要用非对称权重解决残差最小化 [2] 。

分位数回归采用加权残差绝对值之和的方法估计参数，其优点体现在以下几方面：首先，它对模型中的随机扰动项不需做任何分布的假定，这样整个回归模型就具有很强的稳健性；其次，分位数回归本身没有使用一个连接函数来描述因变量的均值和方差的相互关系，因此分位数回归有着比较好的弹性性质；第三，分位数回归由于是对所有分位数进行回归，因此对于数据中出现的异常点具有耐抗性；第四，不同于普通的最小二乘回归，分位数回归对于因变量具有单调变换性；最后，分位数回归估计出来的参数具有在大样本理论下的渐进优良性。

4 计算分位数的C#源程序

using System;

namespace Legalsoft.Truffer
{
    /// <summary>
    /// Object for estimating arbitrary quantile values 
    /// from a continuing stream of data values.
    /// </summary>
    public class IQagent
    {
        public const int nbuf = 1000;
        private int nq { get; set; }
        private int nt { get; set; }
        private int nd { get; set; }
        private double[] pval { get; set; }
        private double[] dbuf;
        private double[] qile { get; set; }
        private double q0 { get; set; }
        private double qm { get; set; }

        public IQagent()
        {
            this.nq = 251;
            this.nt = 0;
            this.nd = 0;
            this.pval = new double[nq];
            this.dbuf = new double[nbuf];
            this.qile = new double[nq];
            this.q0 = 1.0e99;
            this.qm = -1.0e99;

            for (int j = 85; j <= 165; j++)
            {
                pval[j] = (j - 75.0) / 100.0;
            }

            // Set general purpose array of p - values ranging from 1.0e-6 to 1~1.0e-6.
            // You can change this if you want.
            for (int j = 84; j >= 0; j--)
            {
                pval[j] = 0.87191909 * pval[j + 1];
                pval[250 - j] = 1.0 - pval[j];
            }
        }

        /// <summary>
        /// Assimilate a new value from the stream.
        /// </summary>
        /// <param name="datum"></param>
        public void add(double datum)
        {
            dbuf[nd++] = datum;
            if (datum < q0)
            {
                q0 = datum;
            }
            if (datum > qm)
            {
                qm = datum;
            }
            if (nd == nbuf)
            {
                update();
            }
        }

        /// <summary>
        /// Batch update.
        /// This function is called by add or report and should not be called directly
        /// by the user.
        /// </summary>
        public void update()
        {
            int jd = 0;
            int jq = 1;
            double told = 0.0;
            double tnew = 0.0;
            double[] newqile = new double[nq];

            Sorter.sort(dbuf, nd);

            double qold = q0;
            double qnew = q0;
            qile[0] = newqile[0] = q0;
            qile[nq - 1] = newqile[nq - 1] = qm;
            pval[0] = Math.Min(0.5 / (nt + nd), 0.5 * pval[1]);
            pval[nq - 1] = Math.Max(1.0 - 0.5 / (nt + nd), 0.5 * (1.0 + pval[nq - 2]));
            for (int iq = 1; iq < nq - 1; iq++)
            {
                double target = (nt + nd) * pval[iq];
                if (tnew < target)
                {
                    for (; ; )
                    {
                        if (jq < nq && (jd >= nd || qile[jq] < dbuf[jd]))
                        {
                            qnew = qile[jq];
                            tnew = jd + nt * pval[jq++];
                            if (tnew >= target)
                            {
                                break;
                            }
                        }
                        else
                        {
                            qnew = dbuf[jd];
                            tnew = told;
                            if (qile[jq] > qile[jq - 1])
                            {
                                tnew += nt * (pval[jq] - pval[jq - 1]) * (qnew - qold) / (qile[jq] - qile[jq - 1]);
                            }
                            jd++;
                            if (tnew >= target)
                            {
                                break;
                            }
                            told = tnew++;
                            qold = qnew;
                            if (tnew >= target)
                            {
                                break;
                            }
                        }
                        told = tnew;
                        qold = qnew;
                    }
                }
                //if (tnew == told)
                if (Math.Abs(tnew - told) <= float.Epsilon)
                {
                    newqile[iq] = 0.5 * (qold + qnew);
                }
                else
                {
                    newqile[iq] = qold + (qnew - qold) * (target - told) / (tnew - told);
                }
                told = tnew;
                qold = qnew;
            }
            // qile = newqile;
            qile = Globals.CopyFrom(newqile);
            nt += nd;
            nd = 0;
        }

        /// <summary>
        /// Return estimated p-quantile for 
        /// the data seen so far. (E.g., p D 0:5 for median.)
        /// </summary>
        /// <param name="p"></param>
        /// <returns></returns>
        public double report(double p)
        {
            if (nd > 0)
            {
                update();
            }
            int jl = 0;
            int jh = nq - 1;
            int j;
            while (jh - jl > 1)
            {
                j = (jh + jl) >> 1;
                if (p > pval[j])
                {
                    jl = j;
                }
                else
                {
                    jh = j;
                }
            }
            j = jl;
            double q = qile[j] + (qile[j + 1] - qile[j]) * (p - pval[j]) / (pval[j + 1] - pval[j]);
            return Math.Max(qile[0], Math.Min(qile[nq - 1], q));
        }
    }
}