【问题标题】:uncertainty sensitivity analysis with Pomp package Latin hypercube sampling使用 Pomp 包拉丁超立方采样进行不确定性敏感性分析
【发布时间】:2020-04-21 15:38:22
【问题描述】:

一年前,我使用 Pomp 包进行了不确定性分析。目的是量化每个参数的变化对结果变量 (Rhv) 的影响。我将通过拉丁超立方抽样 (LHS) 进行的不确定性分析与通过稳健的偏秩相关系数 (PRCC) 方法进行的敏感性分析相结合。
现在我尝试更新这项工作,它不再工作了。如果可能的话,我希望得到一些帮助。抱歉,帖子太长了。

例子:

##### 1.explain data
library(pomp)
AvgT <- c(19, 30)
b <- 0.0047 * AvgT + 0.0467
Bm <- c(0.37, 0.67)    
Bh <- c(0.17, 0.73)    
eta <- c(0,0.4)
teta <- c(0.1, 0.57)
gammaH <- c(1/8, 1/3)
gammaHl <- c(1/6, 1/4)
gammav <- c(1/10, 1/7)
Ml <- c(0, 200.104)
Hl <- c(0.1352, 30.1182)
c <- c(0.1, 0.3)
muv <- c(1/30,1/3)

##### 2.do the Latin hypercube sampling method to generate 100 000samples
scalars1 <- sobol(vars=list(b = b, Bm = Bm, Bh = Bh,
                            eta = eta, teta = teta, c = c, muv = muv,
                            gammaH = gammaH, gammaHl = gammaHl, gammav = gammav,
                            Ml = Ml, Hl = Hl, AvgT = AvgT), 
                  1e+05)

我现在收到了这条消息,而它在一年前运行良好。我知道 Pomp 包已经更新了。

sobol(vars = list(b = b, Bm = Bm, Bh = Bh, eta = eta, teta = teta, : 找不到函数“sobol”

之后我通常会跑步:

##### 3.create the function
Rhv <- function(b, Bm, Bh, eta, teta, c, gammaH, gammaHl, gammav, muv,
               Ml, Hl, AvgT) {
  sqrt((((b * b * eta * Bm * Bh * ((Ml * c) / (Hl)) * teta) * (1 / gammaH * 1 / muv))
        + ((b * b * Bm * Bh * ((Ml * c) / Hl) * teta) * (1 / gammaHl * 1 / muv))) *
       ((1 / gammav) / (((1 / gammav + 1 / muv)))))
  }
##### 4. do.call allows to call function with names as list 
vals1<- do.call(Rhv, scalars1)

##### 5. bind the columns with scalars results and R0
Both1 <- cbind(scalars1, Rhv=vals1)
head(Both1)
[example data Both1][1]

#####6.PRCC analysis
library('epiR')

然后是一个箱线图来获得我需要的最终数字: example of boxplot

【问题讨论】:

    标签: r


    【解决方案1】:

    sobol 作为独立函数已在几年前从包pomp 中删除,版本为0.65-1。您必须改用函数sobolDesign

    查看在线帮助 (?pomp::sobolDesign) 或 pomp package vignette 部分 design

    【讨论】:

      猜你喜欢
      • 1970-01-01
      • 2017-04-24
      • 1970-01-01
      • 1970-01-01
      • 1970-01-01
      • 2020-02-18
      • 1970-01-01
      • 1970-01-01
      • 1970-01-01
      相关资源
      最近更新 更多