1

Something about how I've written this optimize function is not letting it refer to existing values in my "Input" data set to fill in a new column, "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]]

I get the error

invalid function value in 'optimize'

Here is a brief definition of "KozakTaper" so you have some context for what I'm trying to do. The trick, and why I need to use optimize or something like it, is that I don't want y (which is what KozakTaper returns). I want to know what DHT is if y=10.16, but because the equation can't be rearranged to solve for DHT, I'm using optimize to return the DHT value that minimizes the difference between y and 10.16.

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))}

If you have a different approach for finding DHT for each row of data, I'm open to other suggestions. The optimize strategy works fine for one data point (a tree of given species, diameter, and height), but since I've put in the references to column names, it's hitting the error. All advice on what to modify is appreciated!

-KB

Research Fellow, Appalachian Mountain Club

A simplified example of the input data:

> 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")
thatsneat
  • 25
  • 5
  • It might be easier to troubleshoot this if you provided some sample data in an easily consumed format. Perhaps the output from `dput(head(Input))`. Ref: Refs: https://stackoverflow.com/questions/5963269 – r2evans Oct 29 '18 at 20:52
  • Thanks, @r2evans. I've added a mini data set to the post and am happy to provide more. – thatsneat Oct 30 '18 at 13:51
  • `abs(x)` is non-differentiable. This may cause problems: most solvers expect smooth functions. – Erwin Kalvelagen Nov 04 '18 at 06:24

1 Answers1

0

Your function f returns a vector, not a scalar, as you apply it to the whole data frame. optimize cannot operate on vector-valued functions, thus outputs the message "invalid function value".

Instead, I assume you want to find an optimizing solution for every line of the data frame. Then just do it.

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

(I changed '&' to '&&' in the KozakTaper function to avoid warnings.)

Hans W.
  • 1,799
  • 9
  • 16