R optim 与 Scipy 优化之间的差异：Nelder-Mead答案

【问题标题】：Discrepancies between R optim vs Scipy optimize: Nelder-MeadR optim 与 Scipy 优化之间的差异：Nelder-Mead
【发布时间】：2019-07-25 21:59:12
【问题描述】：

我编写了一个脚本，我认为它应该在 Python 和 R 中产生相同的结果，但它们产生的答案却截然不同。每个尝试通过使用 Nelder-Mead 最小化偏差来将模型拟合到模拟数据。总体而言，R 中的 optim 性能要好得多。难道我做错了什么？ R 和 SciPy 中实现的算法是否不同？

Python 结果：

>>> res = minimize(choiceProbDev, sparams, (stim, dflt, dat, N), method='Nelder-Mead')

 final_simplex: (array([[-0.21483287, -1.        , -0.4645897 , -4.65108495],
       [-0.21483909, -1.        , -0.4645915 , -4.65114839],
       [-0.21485426, -1.        , -0.46457789, -4.65107337],
       [-0.21483727, -1.        , -0.46459331, -4.65115965],
       [-0.21484398, -1.        , -0.46457725, -4.65099805]]), array([107.46037865, 107.46037868, 107.4603787 , 107.46037875,
       107.46037875]))
           fun: 107.4603786452194
       message: 'Optimization terminated successfully.'
          nfev: 349
           nit: 197
        status: 0
       success: True
             x: array([-0.21483287, -1.        , -0.4645897 , -4.65108495])

R 结果：

> res <- optim(sparams, choiceProbDev, stim=stim, dflt=dflt, dat=dat, N=N,
             method="Nelder-Mead")

$par
[1] 0.2641022 1.0000000 0.2086496 3.6688737

$value
[1] 110.4249

$counts
function gradient 
     329       NA 

$convergence
[1] 0

$message
NULL

我检查了我的代码，据我所知，这似乎是由于优化和最小化之间的一些差异，因为我试图最小化的函数（即，choiceProbDev）在每个函数中的操作都是相同的（除了输出，我还检查了函数中每个步骤的等效性）。例如：

Python 选择ProbDev：

>>> choiceProbDev(np.array([0.5, 0.5, 0.5, 3]), stim, dflt, dat, N)
143.31438613033876

R 选择ProbDev：

> choiceProbDev(c(0.5, 0.5, 0.5, 3), stim, dflt, dat, N)
[1] 143.3144

我也尝试过调整每个优化函数的容差水平，但我不完全确定容差参数如何在两者之间匹配。无论哪种方式，到目前为止，我的摆弄并没有使两者达成一致。这是每个的完整代码。

Python：

# load modules
import math
import numpy as np
from scipy.optimize import minimize
from scipy.stats import binom

# initialize values
dflt = 0.5
N = 1

# set the known parameter values for generating data
b = 0.1
w1 = 0.75
w2 = 0.25
t = 7

theta = [b, w1, w2, t]

# generate stimuli
stim = np.array(np.meshgrid(np.arange(0, 1.1, 0.1),
                            np.arange(0, 1.1, 0.1))).T.reshape(-1,2)

# starting values
sparams = [-0.5, -0.5, -0.5, 4]


# generate probability of accepting proposal
def choiceProb(stim, dflt, theta):

    utilProp = theta[0] + theta[1]*stim[:,0] + theta[2]*stim[:,1]  # proposal utility
    utilDflt = theta[1]*dflt + theta[2]*dflt  # default utility
    choiceProb = 1/(1 + np.exp(-1*theta[3]*(utilProp - utilDflt)))  # probability of choosing proposal

    return choiceProb

# calculate deviance
def choiceProbDev(theta, stim, dflt, dat, N):

    # restrict b, w1, w2 weights to between -1 and 1
    if any([x > 1 or x < -1 for x in theta[:-1]]):
        return 10000

    # initialize
    nDat = dat.shape[0]
    dev = np.array([np.nan]*nDat)

    # for each trial, calculate deviance
    p = choiceProb(stim, dflt, theta)
    lk = binom.pmf(dat, N, p)

    for i in range(nDat):
        if math.isclose(lk[i], 0):
            dev[i] = 10000
        else:
            dev[i] = -2*np.log(lk[i])

    return np.sum(dev)


# simulate data
probs = choiceProb(stim, dflt, theta)

# randomly generated data based on the calculated probabilities
# dat = np.random.binomial(1, probs, probs.shape[0])
dat = np.array([0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 1,
       0, 1, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1,
       0, 0, 1, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1,
       0, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1,
       0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1])

# fit model
res = minimize(choiceProbDev, sparams, (stim, dflt, dat, N), method='Nelder-Mead')

R：

library(tidyverse)

# initialize values
dflt <- 0.5
N <- 1

# set the known parameter values for generating data
b <- 0.1
w1 <- 0.75
w2 <- 0.25
t <- 7

theta <- c(b, w1, w2, t)

# generate stimuli
stim <- expand.grid(seq(0, 1, 0.1),
                    seq(0, 1, 0.1)) %>%
  dplyr::arrange(Var1, Var2)

# starting values
sparams <- c(-0.5, -0.5, -0.5, 4)

# generate probability of accepting proposal
choiceProb <- function(stim, dflt, theta){
  utilProp <- theta[1] + theta[2]*stim[,1] + theta[3]*stim[,2]  # proposal utility
  utilDflt <- theta[2]*dflt + theta[3]*dflt  # default utility
  choiceProb <- 1/(1 + exp(-1*theta[4]*(utilProp - utilDflt)))  # probability of choosing proposal
  return(choiceProb)
}

# calculate deviance
choiceProbDev <- function(theta, stim, dflt, dat, N){
  # restrict b, w1, w2 weights to between -1 and 1
  if (any(theta[1:3] > 1 | theta[1:3] < -1)){
    return(10000)
  }

  # initialize
  nDat <- length(dat)
  dev <- rep(NA, nDat)

  # for each trial, calculate deviance
  p <- choiceProb(stim, dflt, theta)
  lk <- dbinom(dat, N, p)

  for (i in 1:nDat){
    if (dplyr::near(lk[i], 0)){
      dev[i] <- 10000
    } else {
      dev[i] <- -2*log(lk[i])
    }
  }
  return(sum(dev))
}

# simulate data
probs <- choiceProb(stim, dflt, theta)

# same data as in python script
dat <- c(0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 1,
         0, 1, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1,
         0, 0, 1, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1,
         0, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1,
         0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
         1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1)

# fit model
res <- optim(sparams, choiceProbDev, stim=stim, dflt=dflt, dat=dat, N=N,
             method="Nelder-Mead")

更新：

在每次迭代打印估计值后，现在在我看来，差异可能源于每种算法采用的“步长”差异。 Scipy 似乎采取了比 optim 更小的步骤（并且在不同的初始方向上）。我还没想好怎么调整。

Python：

>>> res = minimize(choiceProbDev, sparams, (stim, dflt, dat, N), method='Nelder-Mead')
[-0.5 -0.5 -0.5  4. ]
[-0.525 -0.5   -0.5    4.   ]
[-0.5   -0.525 -0.5    4.   ]
[-0.5   -0.5   -0.525  4.   ]
[-0.5 -0.5 -0.5  4.2]
[-0.5125 -0.5125 -0.5125  3.8   ]
...

R：

> res <- optim(sparams, choiceProbDev, stim=stim, dflt=dflt, dat=dat, N=N, method="Nelder-Mead")
[1] -0.5 -0.5 -0.5  4.0
[1] -0.1 -0.5 -0.5  4.0
[1] -0.5 -0.1 -0.5  4.0
[1] -0.5 -0.5 -0.1  4.0
[1] -0.5 -0.5 -0.5  4.4
[1] -0.3 -0.3 -0.3  3.6
...

【问题讨论】：

这是一个非常复杂的用例，而不是简单的optim。我建议打印出数据操作过程的每个步骤，并检查每个数据片段是否在两个脚本中匹配。任何这些步骤之一都可能是问题所在。一张支票：binom.pmf 和 dbinom。
@Parfait -- 谢谢！我已经编辑了我的问题，希望有助于突出和查明问题。我已经编辑了我的代码以更清楚地指定正在使用的库。我逐行浏览了这两个代码，一切似乎都是等效的（例如，我最小化的choiceProbDev 函数在两个实现中产生相同的结果）。
尝试使用 R 的 optim control list 中的 args 来匹配 Python 的 minimize(method='Nelder-Mead') 中的默认值。我认为这是由于默认差异造成的。
Method Nelder-Mead in optim 不是 R 中 Nelder-Mead 的最佳或最准确的实现。您可能想尝试包中的 nmk[b] dfoptim 或 adagio 中的 neldermead[b] 或 pracma 中的自适应版本 anms 等。实现并没有太大的不同，但会导致准确性和效率的显着差异，尤其是在存在多个最小值的情况下。

标签： python r scipy

【解决方案1】：

这并不完全是“优化器差异是什么”的答案，但我想在这里对优化问题进行一些探索。几个要点：

表面是光滑的，因此基于导数的优化器可能会更好地工作（即使没有明确编码的梯度函数，即回归有限差分近似 - 使用梯度函数会更好）
这个表面是对称的，所以它有多个最优值（显然是两个），但它不是高度多峰或粗糙的，所以我认为随机全局优化器不值得麻烦
对于不是太高维或计算成本不太高的优化问题，可视化全局表面以了解正在发生的事情是可行的。
对于有边界的优化，通常最好或者使用显式处理边界的优化器，或将参数的比例更改为不受约束的比例

这是整个表面的图片：

红色轮廓是对数似然等于 (110, 115, 120) 的轮廓（我能得到的最佳拟合是 LL=105.7）。最佳点位于第二列第三行（由 L-BFGS-B 实现）和第五列第四行（真实参数值）。（我还没有检查目标函数以了解对称性来自何处，但我认为这可能很清楚。）Python 的 Nelder-Mead 和 R 的 Nelder-Mead大约同样糟糕。 p>

参数和问题设置

## initialize values
dflt <- 0.5; N <- 1
# set the known parameter values for generating data
b <- 0.1; w1 <- 0.75; w2 <- 0.25; t <- 7
theta <- c(b, w1, w2, t)
# generate stimuli
stim <- expand.grid(seq(0, 1, 0.1), seq(0, 1, 0.1))
# starting values
sparams <- c(-0.5, -0.5, -0.5, 4)
# same data as in python script
dat <- c(0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 1,
         0, 1, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1,
         0, 0, 1, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1,
         0, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1,
         0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
         1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1)

目标函数

注意使用内置函数（plogis()、dbinom(...,log=TRUE)，如果可能的话。

# generate probability of accepting proposal
choiceProb <- function(stim, dflt, theta){
    utilProp <- theta[1] + theta[2]*stim[,1] + theta[3]*stim[,2]  # proposal utility
    utilDflt <- theta[2]*dflt + theta[3]*dflt  # default utility
    choiceProb <- plogis(theta[4]*(utilProp - utilDflt))  # probability of choosing proposal
    return(choiceProb)
}
# calculate deviance
choiceProbDev <- function(theta, stim, dflt, dat, N){
  # restrict b, w1, w2 weights to between -1 and 1
    if (any(theta[1:3] > 1 | theta[1:3] < -1)){
        return(10000)
    }
    ## for each trial, calculate deviance
    p <-  choiceProb(stim, dflt, theta)
    lk <-  dbinom(dat, N, p, log=TRUE)
    return(sum(-2*lk))
}
# simulate data
probs <- choiceProb(stim, dflt, theta)

模型拟合

# fit model
res <- optim(sparams, choiceProbDev, stim=stim, dflt=dflt, dat=dat, N=N,
             method="Nelder-Mead")
## try derivative-based, box-constrained optimizer
res3 <- optim(sparams, choiceProbDev, stim=stim, dflt=dflt, dat=dat, N=N,
              lower=c(-1,-1,-1,-Inf), upper=c(1,1,1,Inf),
             method="L-BFGS-B")

py_coefs <- c(-0.21483287,  -0.4645897 , -1, -4.65108495) ## transposed?
true_coefs <- c(0.1, 0.25, 0.75, 7)  ## transposed?
## start from python coeffs
res2 <- optim(py_coefs, choiceProbDev, stim=stim, dflt=dflt, dat=dat, N=N,
             method="Nelder-Mead")

探索对数似然面

cc <- expand.grid(seq(-1,1,length.out=51),
                  seq(-1,1,length.out=6),
                  seq(-1,1,length.out=6),
                  seq(-8,8,length.out=51))
## utility function for combining parameter values
bfun <- function(x,grid_vars=c("Var2","Var3"),grid_rng=seq(-1,1,length.out=6),
                 type=NULL) {
    if (is.list(x)) {
        v <- c(x$par,x$value)
    } else if (length(x)==4) {
        v <- c(x,NA)
    }
    res <- as.data.frame(rbind(setNames(v,c(paste0("Var",1:4),"z"))))
    for (v in grid_vars)
        res[,v] <- grid_rng[which.min(abs(grid_rng-res[,v]))]
    if (!is.null(type)) res$type <- type
    res
}

resdat <- rbind(bfun(res3,type="R_LBFGSB"),
                bfun(res,type="R_NM"),
                bfun(py_coefs,type="Py_NM"),
                bfun(true_coefs,type="true"))

cc$z <- apply(cc,1,function(x) choiceProbDev(unlist(x), dat=dat, stim=stim, dflt=dflt, N=N))
library(ggplot2)
library(viridisLite)
ggplot(cc,aes(Var1,Var4,fill=z))+
    geom_tile()+
    facet_grid(Var2~Var3,labeller=label_both)+
    scale_fill_viridis_c()+
    scale_x_continuous(expand=c(0,0))+
    scale_y_continuous(expand=c(0,0))+
    theme(panel.spacing=grid::unit(0,"lines"))+
    geom_contour(aes(z=z),colour="red",breaks=seq(105,120,by=5),alpha=0.5)+
    geom_point(data=resdat,aes(colour=type,shape=type))+
    scale_colour_brewer(palette="Set1")

ggsave("liksurf.png",width=8,height=8)

【讨论】：

【解决方案2】：

'Nelder-Mead' 一直是一个有问题的优化方法，它在optim 中的编码不是最新的。我们将尝试 R 包中提供的其他三种实现。

为了避免其他参数，我们将函数fn定义为

fn <- function(theta)
        choiceProbDev(theta, stim=stim, dflt=dflt, dat=dat, N=N)

那么求解器dfoptim::nmk()、adagio::neldermead() 和pracma::anms() 都将返回相同的最小值xmin = 105.7843，但例如在不同的位置

dfoptim::nmk(sparams, fn)
## $par
## [1] 0.1274937 0.6671353 0.1919542 8.1731618
## $value
## [1] 105.7843

这些是真正的局部最小值，而例如 c(-0.21483287,-1.0,-0.4645897,-4.65108495) 处的 Python 解决方案 107.46038 则不是。您的问题数据显然不足以拟合模型。

您可以尝试使用全局优化器，以在特定范围内找到全局最优值。在我看来，所有局部最小值都具有相同的最小值。

【讨论】：

感谢全局优化器的建议。也感谢其他 R 包的有用建议，这些建议可能会更好。如果 optim 已经过时了，那么有趣的是，在这种情况下，它的性能比 scipy.optimize 好得多。 scipy 是不是更过时了？我想我的问题仍然存在：为什么 scipy 的结果与此处列出的所有 R 优化器相比如此不同且差得多？我想了解造成差异的原因，以及是否有办法提高其性能/使其与 R 更接近。
另外，还有另一个 python 模块可以更好地执行 Nelder-Mead 吗？每当我在 python 中搜索优化函数时，似乎绝大多数结果都会导致 scipy.optimize。当然，python 有一个类似的算法，至少可以接近 R 的。
为什么你认为optim 的表现比scipy.minimize 好很多？ optim 至少返回 110.4249，minimize 107.4604，至少使用默认选项。两者都不是真正的局部最小值。 SciPy 源代码只提到了 Nelder 和 Mead 的原始文章，以及 1996 年的概述文章。将自适应 Nelder-Mead 程序（大约 20 LoC）从 R 转换（分别是其 Matlab 版本）并不难) 到 Python。
看看 PyPi 项目nelder-mead，2018 年 10 月的最新版本，这可能是为 Python 社区提供新实现的尝试。使用 Python3 可以使用pip install nelder-mead 安装它。（免责声明：我没有尝试过。）
啊，你是对的...我一直非常关注optim 一直在返回的参数值，这些参数值明显接近于我忽略的原始输入事实上，scipy.optimize 的最小值实际上更低。哎呀。感谢 PyPi 提示。