通过 csv 的行迭代优化函数
Iterating an optimize function through rows of a csv
关于我如何编写这个优化函数的一些事情是不让它引用我的 "Input" 数据集中的现有值来填充新列,"Opt":
Input=read.csv("....csv")
Input$Opt=0
Input$Opt <- optimize(f = function(x) abs(10.16 - KozakTaper(Bark='ob',
SPP=Input$Species,
DHT=x,
DBH=Input$DBH,
HT=Input$Ht,
Planted=0)),
lower=Input$Ht*.25, upper=Input$Ht+1,
maximum = FALSE, tol = .Machine$double.eps^0.25)[[1]]
我收到错误
invalid function value in 'optimize'
这里是 "KozakTaper" 的简要定义,因此您对我正在尝试做的事情有一些了解。诀窍,以及为什么我需要使用优化或类似的东西,是我不想要 y(这就是 KozakTaper returns)。我想知道如果 y=10.16,DHT 是什么,但是因为不能重新排列方程来求解 DHT,所以我使用优化 return 最小化 y 和 10.16 之间差异的 DHT 值。
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))}
如果您有不同的方法来为每行数据查找 DHT,我愿意接受其他建议。优化策略适用于一个数据点(给定物种、直径和高度的树),但由于我已经输入了对列名的引用,所以它遇到了错误。
感谢所有关于修改内容的建议!
-KB
阿巴拉契亚山脉俱乐部研究员
输入数据的简化示例:
> dput(head(Input))
structure(list(Species = structure(c(3L, 3L, 3L, 3L, 8L, 8L), .Label = c("AB",
"BC", "BF", "EH", "PB", "PR", "RM", "RS", "SM", "ST", "WA", "WP",
"YB"), class = "factor"), DBH = c(6.9000001, 8.1000004, 5.8000002,
6.5999999, 9.5, 7.5999999), Ht = c(44, 43, 34, 41, 56, 58)), .Names = c("Species",
"DBH", "Ht"), row.names = c(NA, 6L), class = "data.frame")
您的函数 f
returns 一个矢量,而不是标量,因为您将它应用于整个数据框。 optimize
无法对向量值函数进行运算,因此输出消息 "invalid function value".
相反,我假设您想为数据框的每一行找到一个优化解决方案。那就去做吧。
Bark='ob'; Planted=0
for (i in 1:6) {
SPP=Input$Species[i]; DBH=Input$DBH[i]; HT=Input$Ht[i];
f <- function(x) abs(10.16 - KozakTaper(Bark,SPP,x,DBH,HT,Planted))
o <- optimize(f, lower=Input$Ht[i]*.25, upper=Input$Ht[i]+1,
maximum = FALSE, tol = .Machine$double.eps^0.25)
cat(o$minimum, ' ')
}
## 11.00059 10.75038 8.500041 10.25055 26.2004 26.75473
(我在 KozakTaper 函数中将 '&' 更改为 '&&' 以避免警告。)
关于我如何编写这个优化函数的一些事情是不让它引用我的 "Input" 数据集中的现有值来填充新列,"Opt":
Input=read.csv("....csv")
Input$Opt=0
Input$Opt <- optimize(f = function(x) abs(10.16 - KozakTaper(Bark='ob',
SPP=Input$Species,
DHT=x,
DBH=Input$DBH,
HT=Input$Ht,
Planted=0)),
lower=Input$Ht*.25, upper=Input$Ht+1,
maximum = FALSE, tol = .Machine$double.eps^0.25)[[1]]
我收到错误
invalid function value in 'optimize'
这里是 "KozakTaper" 的简要定义,因此您对我正在尝试做的事情有一些了解。诀窍,以及为什么我需要使用优化或类似的东西,是我不想要 y(这就是 KozakTaper returns)。我想知道如果 y=10.16,DHT 是什么,但是因为不能重新排列方程来求解 DHT,所以我使用优化 return 最小化 y 和 10.16 之间差异的 DHT 值。
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))}
如果您有不同的方法来为每行数据查找 DHT,我愿意接受其他建议。优化策略适用于一个数据点(给定物种、直径和高度的树),但由于我已经输入了对列名的引用,所以它遇到了错误。 感谢所有关于修改内容的建议!
-KB
阿巴拉契亚山脉俱乐部研究员
输入数据的简化示例:
> dput(head(Input))
structure(list(Species = structure(c(3L, 3L, 3L, 3L, 8L, 8L), .Label = c("AB",
"BC", "BF", "EH", "PB", "PR", "RM", "RS", "SM", "ST", "WA", "WP",
"YB"), class = "factor"), DBH = c(6.9000001, 8.1000004, 5.8000002,
6.5999999, 9.5, 7.5999999), Ht = c(44, 43, 34, 41, 56, 58)), .Names = c("Species",
"DBH", "Ht"), row.names = c(NA, 6L), class = "data.frame")
您的函数 f
returns 一个矢量,而不是标量,因为您将它应用于整个数据框。 optimize
无法对向量值函数进行运算,因此输出消息 "invalid function value".
相反,我假设您想为数据框的每一行找到一个优化解决方案。那就去做吧。
Bark='ob'; Planted=0
for (i in 1:6) {
SPP=Input$Species[i]; DBH=Input$DBH[i]; HT=Input$Ht[i];
f <- function(x) abs(10.16 - KozakTaper(Bark,SPP,x,DBH,HT,Planted))
o <- optimize(f, lower=Input$Ht[i]*.25, upper=Input$Ht[i]+1,
maximum = FALSE, tol = .Machine$double.eps^0.25)
cat(o$minimum, ' ')
}
## 11.00059 10.75038 8.500041 10.25055 26.2004 26.75473
(我在 KozakTaper 函数中将 '&' 更改为 '&&' 以避免警告。)