R与Python实现Box-Cox变换所得最优λ值不一致原因咨询

问题背景 #

对同一组数据分别使用三种方式实现Box-Cox变换，得到的最优λ存在明显差异，具体测试代码和结果如下。

测试数据（R格式） #

x <- c(112,118,132,129,121,135,148,148,136,119,104,118,115,126,141,135,125,149,170,170,158,133,114,140,145,150,178,163,172,178,199,199,184,162,146,166,171,180,193,181,183,218,230,242,209,191,172,194,196,196,236,235,229,243,264,272,237,211,180,201,204,188,235,227,234,264,302,293,259,229,203,229,242,233,267,269,270,315,364,347,312,274,237,278,284,277,317,313,318,374,413,405,355,306,271,306,315,301,356,348,355,422,465,467,404,347,305,336,340,318,362,348,363,435,491,505,404,359,310,337,360,342,406,396,420,472,548,559,463,407,362,405,417,391,419,461,472,535,622,606,508,461,390,432)

三种实现方案 #

• 方案A：自定义R脚本

par2 <- 200
par3 <- 100
numlam <- par2 + par3 + 1
n <- length(x)
c <- array(NA,dim=c(numlam))
l <- array(NA,dim=c(numlam))
mx <- -1
mxli <- -999
for (i in 1:numlam)
{
  l[i] <- (i-200-1)/100
  if (l[i] != 0)
  {
    x1 <- (x^l[i] - 1) / l[i]
  } else {
    x1 <- log(x)
  }
  c[i] <- cor(qnorm(ppoints(x), mean=0, sd=1),sort(x1))
  if (mx < c[i])
  {
    mx <- c[i]
    mxli <- l[i]
    x1.best <- x1
  }
}
print(paste("best lambda:", mxli))
hist(x)
hist(x1.best)

• 方案B：R MASS包实现

library(MASS)
bc= boxcox(lm(x~x))
lmbd2 = bc$x[which(bc$y == max(bc$y))]
print(paste("best lambda:", lmbd2))

• 方案C：Python scipy包实现

from scipy.stats import boxcox
import matplotlib.pyplot as plt
x = [112,118,132,129,121,135,148,148,136,119,104,118,115,126,141,135,125,149,170,170,158,133,114,140,145,150,178,163,172,178,199,199,184,162,146,166,171,180,193,181,183,218,230,242,209,191,172,194,196,196,236,235,229,243,264,272,237,211,180,201,204,188,235,227,234,264,302,293,259,229,203,229,242,233,267,269,270,315,364,347,312,274,237,278,284,277,317,313,318,374,413,405,355,306,271,306,315,301,356,348,355,422,465,467,404,347,305,336,340,318,362,348,363,435,491,505,404,359,310,337,360,342,406,396,420,472,548,559,463,407,362,405,417,391,419,461,472,535,622,606,508,461,390,432]
plt.hist(x, bins=10)
y,lmda = boxcox(x)
print(lmda)
plt.hist(y, bins=13)

运行结果 #

三个方案输出的最优λ完全不同：

• 自定义R脚本：0.22
• R MASS包：0.1414
• Python scipy包：0.1480

结果不一致的具体原因 #

三个实现从底层逻辑到计算细节都存在差异，不可能得到完全一致的结果：

• 最优λ的评判标准不统一
  自定义脚本的优化目标是「变换后数据与理论正态分布的相关系数最大」，本质是用正态Q-Q图的线性相关程度选λ，和另外两个包的优化目标完全不同。MASS和scipy虽然都用极大似然估计，但MASS的boxcox是为线性模型场景设计的，计算似然基于线性模型残差，传入的lm(x~x)属于错误的模型设定，会引入额外偏差；scipy是直接对原始独立数据做极大似然估计，不需要依附线性模型框架，两者的目标函数本身就有区别。
• λ的搜索逻辑和精度差异大
  自定义脚本是在[-2, 1]区间按0.01步长做网格遍历；MASS包默认在[-2, 2]区间按0.1步长取粗网格计算似然，直接返回网格点上的最优值，没有做精细化的极值寻优，结果精度很低；scipy使用布伦特一维优化算法直接搜索似然函数的极值点，收敛精度远高于前两种网格搜索方法，结果更接近极大似然的理论最优值。
• 计算细节的处理不一致
  自定义脚本完全没有加入Box-Cox似然估计必须的雅可比校正项，也没有对变换后的数据做中心化处理；MASS和scipy在计算似然时，对方差的无偏校正系数、雅可比项的计算顺序、变换后数据的标准化规则都有细微差异，这些差异累积起来会造成最优λ出现零点零几的偏差。

实际分析中不需要追求λ的精确值，只要λ落在对数似然的95%置信区间内，变换效果就没有本质差异。得到的三个λ对应的变换结果相关度接近1，在稳定方差、正态化的效果上几乎没有区别，完全满足常规统计分析的需求。

内容的提问来源于stack exchange，提问作者David