我试图找到一个逻辑分布函数的不动点,并确定不动点在不同参数值下的变化情况。代码如下所示:
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)
不动点迭代的结构不能找到不稳定的平衡在你的设置,因为它是排斥。换句话说,除非你从不稳定平衡点开始,否则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')
希望这对你有帮助。
我在一个新函数中重新安排了你的代码。
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
我不太明白你用的是哪个发行版;这是我的定点方法的标准代码,我总是使用和更改,如果需要(你必须填写你的函数f(x)在ftn;
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 convergen")
return(NULL)
} else {
return(x.new)
}
}
不知道你到底做错了什么,但我会给你我的代码,这总是为找到固定点工作。下面最后一个函数可以用来计算函数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 convergen")
return(NULL)
} else {
cat("Algorithm convergedn")
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)
}