【问题标题】:Optimize within for loop cannot find function在for循环中优化找不到函数
【发布时间】:2019-03-23 13:52:14
【问题描述】:

我有一个函数,KozakTaper,它返回给定高度 (DHT) 的树干的直径。没有代数方法可以重新排列原始锥度方程以返回给定直径(我的目的是 4 英寸)的 DHT……输入 R! (在 Windows 10 上使用 3.4.3)

我的方法是使用 for 循环迭代 DHT 的可能值(树总高度的 25-100%,HT),然后使用优化选择返回直径最接近 4" 的值。太糟糕了,我得到错误消息Error in f(arg, ...) : could not find function "f"

这是 KozakTaper 的简短定义以及我迄今为止的最佳尝试。

KozakTaper=function(Bark,SPP,DHT,DBH,HT,Planted){
  if(Bark=='ob' & SPP=='AB'){
    a0_tap=1.0693567631
    a1_tap=0.9975021951
    a2_tap=-0.01282775
    b1_tap=0.3921013594
    b2_tap=-1.054622304
    b3_tap=0.7758393514
    b4_tap=4.1034897617
    b5_tap=0.1185960455
    b6_tap=-1.080697381
    b7_tap=0}
  else if(Bark=='ob' & SPP=='RS'){
    a0_tap=0.8758
    a1_tap=0.992
    a2_tap=0.0633
    b1_tap=0.4128
    b2_tap=-0.6877
    b3_tap=0.4413
    b4_tap=1.1818
    b5_tap=0.1131
    b6_tap=-0.4356
    b7_tap=0.1042}
  else{
    a0_tap=1.1263776728
    a1_tap=0.9485083275
    a2_tap=0.0371321602
    b1_tap=0.7662525552
    b2_tap=-0.028147685
    b3_tap=0.2334044323
    b4_tap=4.8569609081
    b5_tap=0.0753180483
    b6_tap=-0.205052535
    b7_tap=0}
  p = 1.3/HT
  z = DHT/HT
  Xi = (1 - z^(1/3))/(1 - p^(1/3))
  Qi = 1 - z^(1/3)
  y = (a0_tap * (DBH^a1_tap) * (HT^a2_tap)) * Xi^(b1_tap * z^4 + b2_tap * (exp(-DBH/HT)) +
                                                    b3_tap * Xi^0.1 + b4_tap * (1/DBH) + b5_tap * HT^Qi + b6_tap * Xi + b7_tap*Planted)
  return(y=round(y,4))}


HT <- .3048*85 #converting from english to metric (sorry, it's forestry)

for (i in c((HT*.25):(HT+1))) {
  d <- KozakTaper(Bark='ob',SPP='RS',DHT=i,DBH=2.54*19,HT=.3048*85,Planted=0)
  frame <- na.omit(d)
  optimize(f=abs(10.16-d), interval=frame, lower=1, upper=90,
           maximum = FALSE,
           tol = .Machine$double.eps^0.25)
}

最终我希望这段代码遍历 csv 并返回 i 以获得最佳 d,这需要重新排列,但我认为我应该首先让它适用于一棵树。
当我打印 d 时,我得到了多个值,因此它正在迭代 i,但它在优化函数中被阻止。 定义frame 是我最近的策略,因为 d 最后返回一个NaN,但它可能不是interval 的最佳输入。我试过interval=c((HT*.25):(HT+1)),在for循环中定义KozakTaper,并在优化之前定义f,但我得到了同样的错误。感谢您对我应该针对的部分(或其他方法)提出建议!

-KB

阿巴拉契亚山脉俱乐部林业研究员。 缅因大学硕士

**使用后续问题进行编辑: 我现在正在尝试为 csv 的每一行“输入”运行此脚本。该行包含 KozakTaper 的值,我用这个来调用它们:

Input=read.csv...
Input$Opt=0
o <- optimize(f = function(x) abs(10.16 - KozakTaper(Bark='ob',
                                                     SPP='Input$Species',
                                                     DHT=x,
                                                     DBH=(2.54*Input$DBH),
                                                     HT=(.3048*Input$Ht),
                                                     Planted=0)),
              lower=Input$Ht*.25, upper=Input$Ht+1,
              maximum = FALSE,  tol = .Machine$double.eps^0.25)
Input$Opt <- o$minimum
Input$Mht <- Input$Opt/.3048. # converting back to English 

Input$Ht 和 Input$DBH 是数字; Input$Species 是因子。 但是,我收到错误invalid function value in 'optimize'。无论我定义“o”还是只运行优化,我都知道。奇怪的是,当我不调用行中的值而是使用答案中的代码时,它告诉我object 'HT' not found。我有一种糟糕的感觉,这是由于我的一些明显/粗心的错误,但我没有通过优化找到有关此错误的帖子。如果您发现我做错了什么,我们将不胜感激!

【问题讨论】:

  • optimize 中的 f= 参数需要一个裸函数:在本例中为 optimize(f=abs, ...)。通过为abs 提供参数,您实际上是在运行abs 函数,然后将其返回值传递给f=
  • 谢谢@divibisan!

标签: r function loops csv mathematical-optimization


【解决方案1】:

我不是优化专家,但我发现三个问题:1) 您对KozakTaper 的调用不会遍历您在循环中指定的范围。 2) KozakTaper 返回单个数字而不是向量。 3)您没有给出优化函数,而是一个表达式。

所以发生的事情是你没有给 optimize 任何东西来迭代。

你只需要这个:

optimize(f = function(x) abs(10.16 - KozakTaper(Bark='ob', 
                                    SPP='RS', 
                                    DHT=x, 
                                    DBH=2.54*19, 
                                    HT=.3048*85, 
                                    Planted=0)),
           lower=HT*.25, upper=HT+1,
           maximum = FALSE,  tol = .Machine$double.eps^0.25)

$minimum
[1] 22.67713 ##Hopefully this is the right answer

$objective
[1] 0

Optimize 现在会将xlower 替换为higher,尽量减少差异

【讨论】:

  • 好的,感谢@emilliman5 解释和编写解决方案!它有效,而且答案非常接近,可以将差异归咎于四舍五入(我引用的是 Excel 电子表格,所以我们会说 R 是正确的)。
  • 嗯实际上它不识别 HT。在运行答案之前,我已经定义了 HT,所以没有意识到它每次都指的是同一个数字。有关更多详细信息,请参阅我在原始帖子底部的编辑。
猜你喜欢
  • 2013-11-15
  • 2011-08-30
  • 1970-01-01
  • 1970-01-01
  • 2017-01-17
  • 2021-07-17
  • 1970-01-01
  • 2015-04-15
相关资源
最近更新 更多