0

4 つの独立変数 (x1、x2、x3、および x4) と 1 つの応答変数 (y、0 から 1 の間の任意の数) を持つデータ セットを分析するより効率的な方法を見つけようとしています。6 つのパラメーター (3 つの勾配 [m1、m2、および m3] と 3 つの切片 [b1、b2、および b3]) を持つモデルにデータを当てはめようとしています。

各ポイントでの y 応答は、次の方程式がゼロに等しい y を解くことによって計算されます。

-x4+(x1/((log(-y/(y-1))-b1)/m1))
          +(x2/((log(-y/(y-1))-b2)/m2))
          +(x3/((log(-y/(y-1))-b3)/m3))

私の知る限り (より良い方法があれば訂正してください)、これを達成する最善の方法は、上記の式の絶対値を で最小化することoptimize()です。以下は、任意のパラメーターと x 値を使用した再現可能な例です。

#model parameters
b1=-8
b2=-10
b3=-15
m1=2
m2=25
m3=50

#independent variables
x1=3.9
x2=.02
x3=.01
x4=1

a=function(y) abs(-x4+(x1/((log(-y/(y-1))-b1)/m1))
      +(x2/((log(-y/(y-1))-b2)/m2))
      +(x3/((log(-y/(y-1))-b3)/m3)))

y=optimize(a,c(0,1))

これらの入力を使用すると、y は ~0.617 と評価されます。

簡単ですが、500 個のデータ ポイントごとにこれを行う必要があります (それぞれに x1、x2、x3、x4 の固有の組み合わせがありますが、パラメータはすべて同じ b1、b2、b3、m1、m2、および m3 です)。私は現在 を使用してvapply()いますが、もっと効率的な方法が必要なようです。ただし、最適化の問題をベクトル化する方法を知りません。

ここで終わってしまえば悪くない（500点全てが1秒以内で評価される）vapply()。optimization()ただし、これは、指定された 1 つのパラメーター セットの y 変数のみを計算します。DEoptim()対数尤度を最大化するためにパラメータ b1、b2、b3、m1、m2、および m3 を最適化しようとすると、問題が発生します(次にDEoptim()、パラメータ微調整の開始パラメータとしてのパラメータを使用しますoptim())。言うまでもなく、パラメータの反復ごとに 500 の最適化問題を評価する必要があるため、これには時間がかかります。奇妙なことに、1 セットのパラメーターを最適化するために手動で vapply を実行すると (つまり、DEoptim の 1 回の反復に相当)、1 秒もかかりませんが、DEoptim()(これも 1 秒もかからないはずです) 約 10 秒かかります。なぜそんなに時間がかかるのかわかりませんがvapply()(何か考えはありますか?)、500 の方程式すべてをより効率的に解く方法があることを願っています。

vapply()とに頼る前に、 ( の代わりにoptimize()) を使用して y の反復計算 (すべての 500 ポイントで同時に) を含むいくつかの代替方法を試しましたが、少し高速でした (反復あたり約 6 秒) が、そうすべきだと思います。反復ごとに 1 秒未満を達成できます (手動で 1 つのインスタンスを実行するときのように)。より良い代替手段が見つからない場合は、おそらく に戻ります。Reduce()optimize()vapply()Reduce()

どんな助けや洞察も大歓迎です。ありがとう！

編集:明確にするために、指定された方程式は、指定されたパラメーターと x 値 (x1、x2、x3、および x4 の一意の組み合わせで構成される合計 500 データ ポイントの合計) の各ポイントで予測された y 値 (方程式がゼロに等しい) を提供します。計算する 500 の一意の y 値があります; y は各インスタンスで唯一不明です)。全体的な目標は、観測された y 値に最適な予測 y 値を取得する最尤法によってパラメーターを最適化することです。m1 m2 m3 および b1 b2 b3 パラメータが変更されているため、パラメータの最適化を繰り返すたびに、500 の y 値を再計算する必要があります。

4

1 に答える 1

3

abs() の代わりに目的関数として正方形を使用することをお勧めします。これは、最初の導関数が連続であるためです。あなたはこれを見ることができます

curve(a,.001,.999)

abs を省略してこの関数を定義することにより、

b=function(y) -x4+(x1/((log(-y/(y-1))-b1)/m1))+
              (x2/((log(-y/(y-1))-b2)/m2))+
              (x3/((log(-y/(y-1))-b3)/m3))

この関数をグラフ化する

curve(b,0.001,.9999)

一般的に言えば、最適化アルゴリズムを使用して連立方程式の解を見つけることは良い考えではありません。これは、最小値 (グローバルおよびローカル) を探すためです。0 の最小値が必要です。

したがって、非線形方程式ソルバーを使用することをお勧めします。それを行うパッケージがありますnleqslv（注：私はそのパッケージの作成者です）。

Vectorize関数は既にベクトル化されているため、ベクトル化の他の方法を使用する必要はありません。

関数を定義します（なしとf同じで、わずかに効率的です）aabs

f <- function(y) {
    tmp <- log(-y/(y-1))
    -x4+(x1/((tmp-b1)/m1))+(x2/((tmp-b2)/m2))+(x3/((tmp-b3)/m3))

}

ヤコビアンを計算する関数を定義します

fjac <- function(x) { h <- 0.00001*x; diag((f(x+h)-f(x))/h) }

の戻り値の各要素fのみが入力ベクトルの対応する要素に依存するため、これは非常に効率的に実行できますy。

パラメータ構成ごとに、データ ベクトルとyソリューションの開始値は次のyように計算できます。

z <- nleqslv(ystart,f,fjac, method="Newton")
y <- z$x

Newtonこの場合、メソッドは機能しないため、メソッドを使用する必要がありますBroyden。この例を試すことができます

K <- 500
x1 <- x1 + c(0,runif(K-1,.1*x1,.3*x1))
x2 <- x2 + c(0,runif(K-1,.01*x2,.03*x2))
x3 <- x3 + c(0,runif(K-1,.01*x3,.03*x3))
x4 <- x4 + c(0,runif(K-1,.01*x4,.03*x4))
nleqslv(rep(.3,K),f,fjac, method="Newton")

私のコンピュータでは、これには約 0.08 秒かかります。

于 2013-08-25T09:36:57.643 に答える