4

私の目標は、分数多項式(FP) の係数を受け入れ、指定された入力数値に対して指定された FP を評価するベクトル化された関数を返す関数を R で記述することです。FP 定義には 2 つの重要なルールがあります。

  • x^0 は log(x) として定義されます
  • べき乗は複数の係数を持つことができます。ここで、べき乗 p の 2 番目の係数は加法項 (x^p*log(x)) に log(x) の係数を追加し、3 番目は log(x)^2 の係数を追加します (x ^p*log(x)^2) など

以下の私の現在のソリューションは、FP 関数を文字列として構築し、文字列を解析して、式を評価する関数を返します。私の質問は、回避するより良い/より高速な方法があるかどうかですeval(parse())-おそらくいくつかのsubstitute()魔法を使用します。

関数は、事前に知られていないが、呼び出されたときに指定されたべき乗あたりの係数の数を処理する必要があります。最終的な FP 評価は頻繁に呼び出されるため、高速である必要があります。

標準の累乗 -2、-1、-0.5、0、0.5、1、2、3 に限定されないことが望ましいでしょう。理想的には、目的の関数は一度に 2 つのステップを実行します: FP 係数を受け入れるだけでなく、数値のベクトルであり、高速でありながら入力の FP 値を返します。

getFP <- function(p_2, p_1, p_0.5, p0, p0.5, p1, p2, p3, ...) {
    p <- as.list(match.call(expand.dots=TRUE)[-1])         # all args
    names(p) <- sub("^p", "", names(p))     # strip "p" from arg names
    names(p) <- sub("_", "-", names(p))     # replace _ by - in arg names

    ## for one power and the i-th coefficient: build string
    getCoefStr <- function(i, pow, coefs) {
        powBT  <- ifelse(as.numeric(pow), paste0("x^(", pow, ")"), "log(x)")
        logFac <- ifelse(i-1,             paste0("*log(x)^", i-1), "")
        paste0("(", coefs[i], ")*", powBT, logFac)
    }

    onePwrStr <- function(pow, p) { # for one power: build string for all coefs
        coefs  <- eval(p[[pow]])
        pwrStr <- sapply(seq(along=coefs), getCoefStr, pow, coefs)
        paste(pwrStr, collapse=" + ")
    }

    allPwrs <- sapply(names(p), onePwrStr, p)  # for each power: build string
    fpExpr  <- parse(text=paste(allPwrs, collapse=" + "))
    function(x) { eval(fpExpr) }
}

例としては-1.5*x^(-1) - 14*log(x) - 13*x^(0.5) + 6*x^0.5*log(x) + 1*x^3、係数 (-1.5, -14, -13, 6, 1) でべき乗 (-1, 0, 0.5, 0.5, 3) を指定したものがあります。

> fp <- getFP(p_1=-1.5, p0=-14, p0.5=c(-13, 6), p3=1)
> fp(1:3)
[1] -13.50000000 -14.95728798   0.01988127
4

1 に答える 1

6

まず、シーケンス内の単一の用語を生成する関数を作成します

one <- function(p, c = 1, repeated = 1) {
  if (p == 0) {
    lhs <- substitute(c * log(x), list(c = c))
  } else {
    lhs <- substitute(c * x ^ p, list(c = c, p = p))
  }

  if (repeated == 1) return(lhs)
  substitute(lhs * log(x) ^ pow, list(lhs = lhs, pow = repeated - 1))
}
one(0)
# 1 * log(x)
one(2)
# 1 * x^2

one(2, 2)
# 2 * x^2

one(2, r = 2)
# 1 * x ^ 2 * log(x)^1
one(2, r = 3)
# 1 * x ^ 2 * log(x)^2

ここで重要なツールは、ここでsubstitute()説明されているものです。

次に、2 つの項を加算する関数を作成します。繰り返しますが、これは代用を使用します:

add_expr_1 <- function(x, y) {
  substitute(x + y, list(x = x, y = y))
}

add_expr_1(one(0, 1), one(2, 1))

これを使用して、任意の数の項を加算する関数を作成できます。

add_expr <- function(x) Reduce(add_expr_1, x)
add_expr(list(one(0, 1), one(1, 1), one(2, 3)))

これらの部品を配置すると、最終的な関数は単純になります。担当者の数を計算し、との組み合わせごとに 1 回Map()呼び出すために使用します。one()powerscoefsreps

fp <- function(powers, coefs) {
  # Determine number of times each power is repeated. This is too
  # clever approach but I think it works
  reps <- ave(powers, powers, FUN = seq_along)

  # Now generate a list of expressions using one
  components <- Map(one, powers, coefs, reps)

  # And combine them together with plus
  add_expr(components)
}

powers <- c(-1, 0, 0.5, 0.5, 3)
coefs <-  c(-1.5, -14, -13, 6, 1)
fp(powers, coefs)
# -1.5 * x^-1 + -14 * log(x) + -13 * x^0.5 + 6 * x^0.5 * log(x)^1 + 
#   1 * x^3
于 2013-11-04T18:14:46.193 に答える