8

ポイントのコレクションで構成されるデータセットがあります。ポイントは、放物線で大まかに囲まれるように平面上に分散されます。ポイントの境界に放物線を合わせる方法を見つけようとしています。

これは私が現在持っているものです：

a = 1
b = 2
c = 3

parabola <- function(x) {
    a * x^2 + b * x + c
}

N = 10000

x <- runif(N, -4, 3)
y <- runif(N, 0, 10)

data <- data.frame(x, y)

data <- subset(data, y >= parabola(x))

plot(data, xlim = c(-5, 5), ylim = c(0, 10), col = "grey")

fr <- function(x) {
    PAR = x[1] * data$x^2 + x[2] * data$x + x[3]
    #
    sum((PAR - data$y)^2 + 100 * plogis(PAR - data$y, scale = 0.00001))
}

par = optim(c(0, 0, 0), fr)$par

a = par[1]
b = par[2]
c = par[3]

curve(parabola, add = TRUE, lty = "dashed")

これにより、サンプルデータセットが作成され、曲線が境界にフィットします。目的関数は、放物線をデータに適合させる「通常の」二乗誤差項と、放物線の下にあるポイントにペナルティを課す2番目のロジスティック項で構成されます。この第2項のパラメーター（100および0.00001）は、試行錯誤によって決定されました。

コードは、点と適合放物線をプロットします。

現在、このシステムは機能しています...しかし、一部の時間のみです。時々それは完全に間違った適合を生成します、そして私はこれらの例ではロジスティック項のパラメータがちょうど不適切であると思います。コードを数回実行して、私が何を意味するかを確認します。

この問題を解決するためのより堅牢な方法が必要だと確信しています。アイデアや提案？

。

4

2 に答える 2

4

私は完全な答えを提供することはできません。私が持っていた唯一のアドホックなアイデアは、最適化アルゴリズムのより良い開始点を提供することでした-最適化しようとする関数の極小値に近づくことを望んでいます。

大まかな最初のバージョンの見積もりはかなり簡単です。b*(x-a)^2+c あなたが見積もることができるようにあなたの放物線を書くならば

a <- data$x[which.min(data$y)]
c <- min(data$y)
 
b1 <- (data$y[which.min(data$x)] - c) / (min(data$x) - a)^2
b2 <- (data$y[which.max(data$x)] - c) / (max(data$x) - a)^2
b <- mean(c(b1, b2))

編集

私の提案と方法「BFGS」を使用して、別の集中的なテストセッションを行いました。次のアプローチでは反例が見つかりませんでした。

seed <- floor(runif(1,1,1000))
set.seed(seed)
a = 1
b = 2
c = 3

parabola <- function(x) {
    b * (x-a)^2 + c
}

N = 10000

x <- runif(N, -4, 3)
y <- runif(N, 0, 10)

data <- data.frame(x, y)

data <- subset(data, y >= parabola(x))

plot(data, xlim = c(-5, 5), ylim = c(0, 10), col = "grey")

fr <- function(x) {
    PAR = x[2] * (data$x - x[1])^2 + x[3]
    #
    sum((PAR - data$y)^2 + 100 * plogis(PAR - data$y, scale = 0.00001))
}

a <- data$x[which.min(data$y)]
c <- min(data$y)

b1 <- (data$y[which.min(data$x)] - c) / (min(data$x) - a)^2
b2 <- (data$y[which.max(data$x)] - c) / (max(data$x) - a)^2
b <- mean(c(b1, b2))

par = optim(c(a, b, c), fr, method="BFGS")$par

a = par[1]
b = par[2]
c = par[3]

curve(parabola, add = TRUE, lty = "dashed")

ただし、正しい収束は保証されません。私は約50のケースを試しましたが、すべてうまくいきました。結果はレビューされていますか、それとも自動化されて正しく機能する必要がありますか？

編集2

目的関数をより信頼性の高いものに更新する方法について、いくつか考えました。今のところ、完全な解決策を考え出す時間はありませんが、おそらくこの考えがあなたを助けるかもしれません：

内に日付がありますrange(data$x)。ここで、このデータの下限にできるだけ適合する放物線を見つけたいと思います。つまり、最大化する値a、b、cを見つけます。

\int_{\range(x)} ax^2 + bx+c dx

（不器用なLaTeXを許してください-数式を書く方が良い場合もあります）。

これで、放物線の下のペナルティポイントは、次のようなペナルティ関数を使用して実行できます。

\lambda (ax_i^2+bx_i+c - y_i)^2 if below parabola, 0 otherwise

区間からその関数を引くと、適切で滑らかな目的関数が得られます。関数を可能な限り単純化することは、データポイントの中央を通る線を近似しようとする最小二乗アプローチを使用するよりも優れたモデルのようです。

ただし、適切なラムダを選択する必要があります。しかし、それは典型的なことです。2つの異なる目的（データの適合、放物線の最大化）の間で妥協する必要があります。どちらがより重要であるかはあなたによって提出されなければなりません。

于 2012-11-23T07:40:31.057 に答える
0

さらに、彼の非常に有益な提案と私の素朴なアイデアの修正をしてくれたthiloに感謝します。放物線の下の領域と適切なペナルティ関数を使用して、thiloの提案に基づいて、以下の解決策が機能するようです。また、Nが小さいほどパフォーマンスが向上するため、L-BFGS-B最適化に変更しました。

parabola.objective <- function(p) {
    d = p[2] * (data$x - p[1])^2 + p[3] - data$y
    #
    area <- function(x) {
        p[2] / 3 * (x - p[1])^3 + p[3] * x
    }
    #
    sum(- area(max(data$x)) + area(min(data$x)) + 100 * ifelse(d > 0, d^2, 0))
}

A <- data$x[which.min(data$y)]
C <- min(data$y)

B1 <- (data$y[which.min(data$x)] - C) / (min(data$x) - A)^2
B2 <- (data$y[which.max(data$x)] - C) / (max(data$x) - A)^2
B <- mean(c(B1, B2))

# the key to getting this working with a small number of points is the
# optimisation method: BFGS works well with around 300 points or more
# but L-BFGS-B seems to perform better down to around 100 points.
#
O = optim(c(A, B, C), parabola.objective, method="L-BFGS-B")

par = O$par

A = par[1]
B = par[2]
C = par[3]

curve(parabola, add = TRUE, lty = "dashed")
于 2012-11-26T07:03:57.750 に答える