【问题标题】:Non linear Constrained Optimisation in RR中的非线性约束优化
【发布时间】:2019-03-28 22:39:23
【问题描述】:

我正在尝试预测特定疾病的存活率。我掌握的唯一信息是诊断后 1 年、3 年、5 年和 10 年的存活率。

例如:

S,即成活率

S<-c(81,78,72,65)

x,诊断后的时间

x<-c(1,3,5,10)

我正在尝试测试一些可以让我估算 20 年后存活率的函数。

我的一个函数被定义为

f(x)= exp(ax^b),a 和 b 未知,但需要为正。我非常友好地使用了 fmarm 提供的代码,但用另一个函数进行了测试。

f(x) = (1 + (x/a)^b)^-1

但是,我得到了非常奇怪的值,都低于 1,我似乎不知道为什么。我错过了什么吗?

S<-c(81,78,72,65)
x<-c(1,3,5,10)

f<-function(ab)
{
  a <- ab[1]
  b <- ab[2]
  return(sum((((1+(x/a)**b)**-1)-S)**2))
}

minim <- nlm(f,p=c(1,1))

ab <- minim$estimate

a_opt <- ab[1]
b_opt <- ab[2]

prediction_exp <- function(x){
  return((1+(x/a_opt)**b_opt)**-1)
}
prediction_exp(20)

plot(prediction_exp(1:20), type="l", col="blue", xlab="Nb d'années après diagnostic", ylab="survie nette en %")
lines(x,S,col="black")

P.S:我发现了我的错误。 S 向量需要小于 1,并且函数应该是 x*a 而不是 (x/a)。再次感谢 fmarm 帮助我!

【问题讨论】:

    标签: r optimization sum


    【解决方案1】:

    在您的情况下,S 和 x 是固定的,您希望找到使 sum(i=1 to 4) exp(a*x[i]**b)-S[i])** 最小化的 a 和 b 2

    你可以创建一个函数

    f <- function(ab){
      a <- ab[1]
      b <- ab[2]
      return(sum((exp(a*x**b)-S)**2))
    }  
    

    ab 是一个长度为 2 的向量,其中 a 在第一位,b 在第二位

    要最小化这个功能,你可以使用nlm

    minim <- nlm(f,p=c(0,0))
    

    你必须给 p : ab 的起始参数。因为我不知道有什么好我只是把 a=0 和 b=0 结果有一个估计组件,可为您提供算法找到的最佳参数

    ab <- minim$estimate
    

    然后你可以从ab中提取a和b

    a_opt <- ab[1]
    b_opt <- ab[2]
    

    你可以创建你的预测函数

    prediction_exp <- function(x){
      return(exp(a_opt*x**b_opt))
    }
    prediction_exp(20)
    

    预计20年后的存活率约为63%

    【讨论】:

    • 非常感谢!
    • SO 设置建议您不要使用 cmets 来表示感谢,而是提供您已正确使用的复选标记系统。就我个人而言,我预计存活率将在 [0-1] 范围内
    【解决方案2】:

    这只是接受的答案中提供的代码,具有适当的生存结果(限制在区间 [0-1] 内,并得到更正的结果:

    S<-c(81,78,72,65)/100
    x<-c(1,3,5,10)
    
    f<-function(ab)
    {
        a <- ab[1]
        b <- ab[2]
        return(sum((((1+(x*a)**b)**-1)-S)**2))
    }
    
    minim <- nlm(f,p=c(1,1))
    
    ab <- minim$estimate
    
    a_opt <- ab[1]
    b_opt <- ab[2]
    
    prediction_exp <- function(x){
        return((1+(x*a_opt)**b_opt)**-1)
    }
    prediction_exp(20)
    [1] 0.5975635
    
    png(); plot(prediction_exp(1:20), type="l", col="blue", xlab="Nb d'années après diagnostic", ylab="survie nette en %")
    lines(x,S,col="black") ; dev.off()
    

    【讨论】:

      猜你喜欢
      • 2014-03-13
      • 2013-02-15
      • 1970-01-01
      • 1970-01-01
      • 1970-01-01
      • 1970-01-01
      • 1970-01-01
      • 1970-01-01
      • 1970-01-01
      相关资源
      最近更新 更多