寻找函数的不动点答案

【问题标题】：Finding the fixed points of a function寻找函数的不动点
【发布时间】：2013-03-06 20:17:57
【问题描述】：

我试图找到一个逻辑分布函数的不动点，并确定不同参数值的不动点如何变化。代码如下：

nfxp.reps <- 0
err <- 10
p <- seq(0, 1, by = 0.0001)
pold <- p
gamma <- 6
k <- 3
while (err > 1E-12) {
  nfxp.reps <- nfxp.reps + 1
  if (nfxp.reps > 999) { 
    stop("The number of NFXP reps needs to be increased. \n")
  } 
  pnew <- plogis(-k + gamma * pold) 
  err <- max(abs(pnew - pold))
  pold <- pnew
}

上面的代码在上面的参数选择中效果很好：gamma 和 k - 找到 3 个固定点，2 个稳定点和 1 个不稳定点（其中 p=0.5）。但是，如果我不按比例更改上述参数，其中中间固定点高于或低于 0.5，例如：

gamma<-7
k<-3

循环无法定位中间不动点 p=0.3225（如果 gamma=7，k=3）

【问题讨论】：

标签： r iteration fixed-point-iteration

【解决方案1】：

构造定点迭代无法在您的设置中找到不稳定的平衡，因为它是排斥的。换句话说，除非您从不稳定的平衡开始，否则 nfxp 算法将始终远离它。

另一种方法是使用根求解方法。当然，不能保证会找到所有固定点。这是一个简单的例子：

library(rootSolve) # for the uniroot.all function
pfind<-function(k=3,gamma=7) 
{
pdiff <-function(p0) p0-plogis(-k + gamma * p0) 
uniroot.all(p.diff,c(0,1))
}
> fps= pfind()
> fps
[1] 0.08036917 0.32257992 0.97925817

我们可以验证这一点：

pseq =seq(0,1,length=100)
plot(x=pseq ,y= plogis(-k + gamma *pseq),type= 'l')
abline(0,1,col='grey')
points(matrix(rep(fps,each=2), ncol=2, byrow=TRUE),pch=19,col='red')

希望这会有所帮助。

【讨论】：

【解决方案2】：

我在一个新函数中重新排列你的代码。

p.fixed <- function(p0,tol = 1E-9,max.iter = 100,k=3,gamma=7,verbose=F){
  pold <- p0
  pnew <-  plogis(-k + gamma * pold) 
  iter <- 1
    while ((abs(pnew - pold) > tol) && (iter < max.iter)){
      pold <- pnew
      pnew <- plogis(-k + gamma * pold) 
      iter <- iter + 1
      if(verbose)
         cat("At iteration", iter, "value of p is:", pnew, "\n")
    }
    if (abs(pnew - pold) > tol) {
      cat("Algorithm failed to converge")
      return(NULL)
    }
    else {
      cat("Algorithm converged, in :" ,iter,"iterations \n")
      return(pnew)
    }
}

一些测试：

p.fixed(0.2,k=3,gamma=7)
Algorithm converged, in : 30 iterations 
[1] 0.08035782
> p.fixed(0.2,k=5,gamma=5)
Algorithm converged, in : 7 iterations 
[1] 0.006927088
> p.fixed(0.2,k=5,gamma=5,verbose=T)
At iteration 2 value of p is: 0.007318032 
At iteration 3 value of p is: 0.006940548 
At iteration 4 value of p is: 0.006927551 
At iteration 5 value of p is: 0.006927104 
At iteration 6 value of p is: 0.006927089 
At iteration 7 value of p is: 0.006927088 
Algorithm converged, in : 7 iterations 
[1] 0.006927088

【讨论】：

感谢以上。但是，如您所见，上面的代码也无法在任何参数规范下找到不稳定的平衡。在 k=3 和 gamma=7 下，固定点应该是 (0.08035, 0.3225, 097926) 它似乎总是错过中间的。
@user1682980 你怎么知道结果？我该如何检查？
只是通过图形测试，最终代码需要在结构中添加大量协变量，所以这只是一个简单的循环开始，这样我就可以确定所有固定点都找到了, 当某些参数值有多个固定点时。

【解决方案3】：

我不太明白您使用的是哪个发行版；这是我经常使用的定点方法的标准代码，如果需要，我会更改（您必须在 ftn 中填写函数 f(x)；

fixed_point <- function(x0, eps = 1e-6, max_iter = 100){
  x.old <- x0
  x.new <- ftn(x.old)
  iter <- 1
  while((abs(x.new-x.old) > eps) && (iter < max_iter){
    x.old <- x.new
    x.new <- ftn(x.old)
    iter <- iter + 1
  }
 if (abs(x.new-x.old) > eps){
  cat("failed to converge\n")
  return(NULL)
 } else {
  return(x.new)
 }
}

【讨论】：

【解决方案4】：

不确定你到底做错了什么，但我会给你我的代码，它总是可以找到固定点。下面最后一个函数可以用来计算函数g，定义为g(x) = c*ftn(x) + x。

fixpt_own <- function(x0, tol = 1e-6, max.iter = 100) {
  xold <- x0
  xnew <- ftn_g(xold)
  iter <- 1
  cat("At iteration 1 value of x is:", xnew, "\n")
  while ((abs(xnew-xold) > tol) && (iter < max.iter)) {
    xold <- xnew;
    xnew <- ftn_g(xold);
    iter <- iter + 1
    cat("At iteration", iter, "value of x is:", xnew, "\n")
  }
  if (abs(xnew-xold) > tol) {
    cat("Algorithm failed to converge\n")
    return(NULL)
  } else {
    cat("Algorithm converged\n")
    return(xnew)
  }
}
fixpt_own(3,1e-6,150)


ftn_g <- function(x){
  c <- 4;
  g <- c*(((1+x)/x - log(x))/(1+x)^2) + x;
  return(g)
}

【讨论】：