Pearson correlation coefficient 的简单理解1. 背景最近LZ正在恶补一些nlp相关的知识其中在判断句子文本相似性的时候使用了Pearson correlation和Spearman’s rank correlation这两个统计学相关系数还是做个笔记吧不然时间久了又会忘记了。在统计学中皮尔逊相关系数Pearson correlation coefficient, PCC,有的时候也指Pearson’s r, 皮尔逊剂矩相关系数 Pearson product-moment correlation coefficient, PPMCC或者双变量相关bivariate correlation,它主要是用来度量两组数据之间线性相关性的。具体实现是两个变量的协方差除以他们的标准差的乘积因此他本质上是对协方差的归一化使得结果始终具有介于-1到1之间的值。与协方差本身一样这个度量只能反映变量的线性相关性而忽略了许多其他类型的关系或相关性。举个简单的例子通常家长都会期望初中高中时期孩子的年龄和身高具有相助大于0且小于1的PCC通常情况下1表示不切实际的完美相关性2. 定义两个变量之间的PCC定义为两个变量的协方差除以他们标准差的乘积ρX,Ycov(X,Y)σXσYE[(X−μX)(Y−μY)]σXσY\rho_{X,Y}\frac {cov(X,Y)}{\sigma_X\sigma_Y}\frac{E[(X-\mu_X)(Y-\mu_Y)]}{\sigma_X\sigma_Y}ρX,YσXσYcov(X,Y)σXσYE[(X−μX)(Y−μY)]上式定义了总体相关系数常用希腊小写字母ρ\rhoρ作为代表符号估算样本的协方差和标准差可以得到样本相关系数PCC,常用英文小写字母rrr表示r∑i1n(Xi−X‾)(Yi−Y‾)∑i1n(Xi−X‾)2∑i1n(Yi−Y‾)2r \frac{\sum_{i1}^{n}{(X_i-\overline{X})(Y_i-\overline{Y})}}{\sqrt{\sum_{i1}^{n}{(X_i-\overline{X})^2}}\sqrt{\sum_{i1}^{n}{(Y_i-\overline{Y})^2}}}r∑i1n(Xi−X)2∑i1n(Yi−Y)2∑i1n(Xi−X)(Yi−Y)r也可以由(Xi,Yi)(X_i, Y_i)(Xi,Yi)样本点的标准分数均值估计得到与上式等价的表达式r1n−1∑i1n(Xi−X‾σX)(Yi−Y‾σY)r\frac{1}{n-1}\sum_{i1}^{n}{(\frac{X_i-\overline{X}}{\sigma_X})(\frac{Y_i-\overline{Y}}{\sigma_Y})}rn−11∑i1n(σXXi−X)(σYYi−Y)其中Xi−X‾σX\frac{X_i-\overline{X}}{\sigma_X}σXXi−X,X‾\overline{X}X,σX\sigma_XσX分别是XiX_iXi样本的标准分数、样本平均值和样本标准差。3. 数学特性数据总体和样本的PCC的绝对值是介于0和1之间如果样本数据点准确的落在直线上计算样本PCC的情况或者双变量分布完全在直线上计算总体PCC的情况则相关系数等于1或者等于-1,PCC是对称的也就是说corr(X,Y)corr(Y,X)corr(X,Y)corr(Y,X)corr(X,Y)corr(Y,X)PCC有一个重要的数学特性是因两个变量的位置和尺度的变化并不会引起该系数的改变也就是说如果我们对X进行abXabXabX和把Y进行cdYcdYcdY进行平移和尺度缩放其中a,b,c,d都是大于0的常数这种操作并不会改变PCC的值4. 代码实现分别使用numpy自己实现对应函数调用scipy中的函数两者结果是一致的自己实现的函数比较简单并没有考虑很多异常情况有兴趣的小伙伴可以去自己看下scipy中的函数实现。fromscipyimportstatsimportnumpyasnpdefComputePCC(x,y):dtypetype(1.0x[0]y[0])xmeanx.mean(dtypedtype)ymeany.mean(dtypedtype)xmx.astype(dtype)-xmean ymy.astype(dtype)-ymean normxmnp.linalg.norm(xm)normymnp.linalg.norm(ym)rnp.dot(xm/normxm,ym/normym)returnr anp.array([0,0,0,1,1,1,1])bnp.arange(7)print(stats.pearsonr(a,b)[0])print(ComputePCC(a,b))defpearsonr(x,y):r Pearson correlation coefficient and p-value for testing non-correlation. The Pearson correlation coefficient [1]_ measures the linear relationship between two datasets. The calculation of the p-value relies on the assumption that each dataset is normally distributed. (See Kowalski [3]_ for a discussion of the effects of non-normality of the input on the distribution of the correlation coefficient.) Like other correlation coefficients, this one varies between -1 and 1 with 0 implying no correlation. Correlations of -1 or 1 imply an exact linear relationship. Positive correlations imply that as x increases, so does y. Negative correlations imply that as x increases, y decreases. The p-value roughly indicates the probability of an uncorrelated system producing datasets that have a Pearson correlation at least as extreme as the one computed from these datasets. Parameters ---------- x : (N,) array_like Input array. y : (N,) array_like Input array. Returns ------- r : float Pearsons correlation coefficient. p-value : float Two-tailed p-value. Warns ----- PearsonRConstantInputWarning Raised if an input is a constant array. The correlation coefficient is not defined in this case, so np.nan is returned. PearsonRNearConstantInputWarning Raised if an input is nearly constant. The array x is considered nearly constant if norm(x - mean(x)) 1e-13 * abs(mean(x)). Numerical errors in the calculation x - mean(x) in this case might result in an inaccurate calculation of r. See Also -------- spearmanr : Spearman rank-order correlation coefficient. kendalltau : Kendalls tau, a correlation measure for ordinal data. Notes ----- The correlation coefficient is calculated as follows: .. math:: r \frac{\sum (x - m_x) (y - m_y)} {\sqrt{\sum (x - m_x)^2 \sum (y - m_y)^2}} where :math:m_x is the mean of the vector :math:x and :math:m_y is the mean of the vector :math:y. Under the assumption that :math:x and :math:m_y are drawn from independent normal distributions (so the population correlation coefficient is 0), the probability density function of the sample correlation coefficient :math:r is ([1]_, [2]_): .. math:: f(r) \frac{{(1-r^2)}^{n/2-2}}{\mathrm{B}(\frac{1}{2},\frac{n}{2}-1)} where n is the number of samples, and B is the beta function. This is sometimes referred to as the exact distribution of r. This is the distribution that is used in pearsonr to compute the p-value. The distribution is a beta distribution on the interval [-1, 1], with equal shape parameters a b n/2 - 1. In terms of SciPys implementation of the beta distribution, the distribution of r is:: dist scipy.stats.beta(n/2 - 1, n/2 - 1, loc-1, scale2) The p-value returned by pearsonr is a two-sided p-value. For a given sample with correlation coefficient r, the p-value is the probability that abs(r) of a random sample x and y drawn from the population with zero correlation would be greater than or equal to abs(r). In terms of the object dist shown above, the p-value for a given r and length n can be computed as:: p 2*dist.cdf(-abs(r)) When n is 2, the above continuous distribution is not well-defined. One can interpret the limit of the beta distribution as the shape parameters a and b approach a b 0 as a discrete distribution with equal probability masses at r 1 and r -1. More directly, one can observe that, given the data x [x1, x2] and y [y1, y2], and assuming x1 ! x2 and y1 ! y2, the only possible values for r are 1 and -1. Because abs(r) for any sample x and y with length 2 will be 1, the two-sided p-value for a sample of length 2 is always 1. References ---------- .. [1] Pearson correlation coefficient, Wikipedia, https://en.wikipedia.org/wiki/Pearson_correlation_coefficient .. [2] Student, Probable error of a correlation coefficient, Biometrika, Volume 6, Issue 2-3, 1 September 1908, pp. 302-310. .. [3] C. J. Kowalski, On the Effects of Non-Normality on the Distribution of the Sample Product-Moment Correlation Coefficient Journal of the Royal Statistical Society. Series C (Applied Statistics), Vol. 21, No. 1 (1972), pp. 1-12. Examples -------- from scipy import stats a np.array([0, 0, 0, 1, 1, 1, 1]) b np.arange(7) stats.pearsonr(a, b) (0.8660254037844386, 0.011724811003954649) stats.pearsonr([1, 2, 3, 4, 5], [10, 9, 2.5, 6, 4]) (-0.7426106572325057, 0.1505558088534455) nlen(x)ifn!len(y):raiseValueError(x and y must have the same length.)ifn2:raiseValueError(x and y must have length at least 2.)xnp.asarray(x)ynp.asarray(y)# If an input is constant, the correlation coefficient is not defined.if(xx[0]).all()or(yy[0]).all():warnings.warn(PearsonRConstantInputWarning())returnnp.nan,np.nan# dtype is the data type for the calculations. This expression ensures# that the data type is at least 64 bit floating point. It might have# more precision if the input is, for example, np.longdouble.dtypetype(1.0x[0]y[0])ifn2:returndtype(np.sign(x[1]-x[0])*np.sign(y[1]-y[0])),1.0xmeanx.mean(dtypedtype)ymeany.mean(dtypedtype)# By using astype(dtype), we ensure that the intermediate calculations# use at least 64 bit floating point.xmx.astype(dtype)-xmean ymy.astype(dtype)-ymean# Unlike np.linalg.norm or the expression sqrt((xm*xm).sum()),# scipy.linalg.norm(xm) does not overflow if xm is, for example,# [-5e210, 5e210, 3e200, -3e200]normxmlinalg.norm(xm)normymlinalg.norm(ym)threshold1e-13ifnormxmthreshold*abs(xmean)ornormymthreshold*abs(ymean):# If all the values in x (likewise y) are very close to the mean,# the loss of precision that occurs in the subtraction xm x - xmean# might result in large errors in r.warnings.warn(PearsonRNearConstantInputWarning())rnp.dot(xm/normxm,ym/normym)# Presumably, if abs(r) 1, then it is only some small artifact of# floating point arithmetic.rmax(min(r,1.0),-1.0)# As explained in the docstring, the p-value can be computed as# p 2*dist.cdf(-abs(r))# where dist is the beta distribution on [-1, 1] with shape parameters# a b n/2 - 1. special.btdtr is the CDF for the beta distribution# on [0, 1]. To use it, we make the transformation x (r 1)/2; the# shape parameters do not change. Then -abs(r) used in cdf(-abs(r))# becomes x (-abs(r) 1)/2 0.5*(1 - abs(r)). (r is cast to float64# to avoid a TypeError raised by btdtr when r is higher precision.)abn/2-1prob2*special.btdtr(ab,ab,0.5*(1-abs(np.float64(r))))returnr,probps 理解还是比较浅显的主要参考的是wiki上的介绍公式中没有使用加粗的形式表示向量和矩阵说明一下哈。
