2016-04-21 12 views
1

私は生物学的集団のダイナミクスをシミュレートする関数が行列内に含まれている数学的モデルを使って作業しています。これは通常、関数を呼び出して行列を更新し、行列の乗算を実行してシステムのダイナミクスを次の時間ステップに投影することを必要とします。私は、可能な場合、マトリックス内に直接関数を埋め込み、明示的に行列要素を更新することをスキップすることができますか?In R関数を含む行列に対して行列演算を行うことは可能ですか?

モデルは次のようになります。 システムの状態が行列X1であり、例えば

X <- c(0,10) 

システムが変わる要素NAが関数である

A <- matrix(data = c(0, NA,0.75,0.75),nrow =2,byrow = T) 
    [,1] [,2] 
[1,] 0.00 NA 
[2,] 0.75 0.75 

をマトリクスに応じてX1の要素2の、このように

f.A.12 <- function(X.2){1 + (1 - X.2/1000)} 

E本: 更新マリックスA更新されたマトリックスとX

A[1,2] <- f.A.12(X[2]) 

反復モデルの状態に基づいて:

A%*%X 

    [,1] 
[1,] 19.9 
[2,] 7.5 

これは、各乗算

本物の前にAを更新し、1000倍を繰り返しますモデルは、複数の関数を含むはるかに大きな行列を使用します。 は私がそのような

A.with.fx <- matrix(data = c(0, f.A.12, 
          0.75, 0.75), nrow =2,byrow = T) 

として、マトリックス内で直接機能を埋め込むためにそして定期的な行列演算を実行することを可能にするRパッケージは、明示的に割り当てることがなくても、例えば

A.with.fx%*%X 

、あります各反復でXの関数であるAの値? これは、独自の必要なルックアップを行う変更された%*%操作である関数が必要になると思います。

+2

すべての要素がなければなりません同じ種類。つまり、関数の行列を持つことはできますが、数値型の行列に関数を埋め込むことはできません。 –

答えて

2

私はあなたが非常にきれいでなく、標準的ではない方法であなたの問題を打破しようとしているかもしれないと思います。スカラ関数を含む行列について考えるのではなく、行列を返す関数について考えてみてください。

あなたは自明行うことができます。

update <- function(X) matrix(c(0,.75, f.A.12(X[2]), .75), 2) 

と表記のいくつかの乱用と:

`%**%` <- function(f, X) f(X) %*% X 
X <- update %**% X 

%*%は一般的なS4であるので、それはそれをオーバーロードするために、痛みのビットだと。

また、少しの代数であなたはXの定数や関数にAを分割し、複雑さをそのように処理することができますので、アップデートを書き換えることができ

A = [ [ 0, 2 - x[2]/1000], 
     [.75, .75]   ] 
    = [[0, 2],[.75, .75]] + [[0, - x[2]/1000],[0, 0 ] ] 
    = c + g(x) 

X = c * X + g(X) * X    

gは行列値関数であり、より遅いがよりクリーナーをもたらす。

g <- function(X) matrix(c(0,0, -x[2]/1000), 0), 2,2) 
C <- matrix(c(0,.75,2,.75),2) 
X <- c(0,10) 

> (X <- C %*% X + g(X) %*% X) 
    [,1] 
[1,] 19.9 
[2,] 7.5 

私たちは一つの関数にg(X) %*% Xを構成することにより、少し良く操作を行うことができます与え

f <- function(x) c(-x[2]^2/1000, 0 ) 

> X <- c(0,10) 
> (X <- C %*% X + f(X)) 
    [,1] 
[1,] 19.9 
[2,] 7.5 

はここでカップルのベンチマークです:マトリックス中

microbenchmark(two_step={A <- update(X); A %*% X}, 
        split=(C %*% X + g(X) %*% X), 
       composed= C %*% X + f(X) ) 
Unit: microseconds 
    expr min  lq mean median  uq max neval 
two_step 3.608 4.0900 4.85369 4.3730 4.9290 13.089 100 
    split 5.587 6.0400 7.72511 6.4900 7.3785 53.047 100 
composed 2.266 2.4835 2.78195 2.6815 2.9745 6.697 100 
関連する問題