ループを回避するベクトル化関数

Question

非常に長く実行されているため、コードをスピードアップしようとしています。私はすでに問題がどこにあるのかを知っていました。次の例を考える：

x<-c((2+2i),(3+1i),(4+1i),(5+3i),(6+2i),(7+2i)) 
P<-matrix(c(2,0,0,3),nrow=2) 
out<-sum(c(0.5,0.5)%*%mtx.exp(P%*%(matrix(c(x,0,0,x),nrow=2)),5)) 

Iが複素数値を有するベクトルxを有し、ベクトルは12^11のエントリを持っており、私は3行目の和を計算します。 （関数はパッケージのBiodemにありますので、関数mtx.expが必要です。％^％関数は複雑な引数をサポートしていません）。

私はエラーを取得

sum(c(0.5,0.5)%*%mtx.exp(P%*%(matrix(c(x,0,0,x),nrow=2)),5)) 

を試してみてください：「ポット％の誤差が*％ポット：非適合引数を」だから私の解決策は、ループを使用することでした：

tmp<-NULL 
for (i in 1:length(x)){ 
    tmp[length(tmp)+1]<-sum(c(0.5,0.5)%*%mtx.exp(P%*%matrix(c(x[i],0,0,x[i]),nrow=2),5)) 
} 

しかし、これは非常に時間がかかります。コードをスピードアップする方法はありますか？私もサプリを試みましたが、それはループと同じくらい長くかかります。

私はこの関数をapproximatly 500回実行しなければならないので、私が助けてくれることを願っています。これは最初に3時間以上かけて試しました。これは非常に満足していない..です

は、コードは、事前に割り当て、あなたのベクトルでスピードアップ

tmp <- rep(NA,length(x)) 

が、私は本当にあなたがしようとしているのか理解していないことができ、非常に

Answer 1

をありがとう計算：最初の例では 、
あなたは（^を行うことができます）対角行列 の電源を取っている、第二で 、非正方行列の電源を取るしようとしています。

次はあなたの計算に相当すると思わ：

sum(P^5/2) * x^5 

EDIT

Pは対角線とCスカラーではないではない場合、 私はmtx.exp(P %*% C, 5)のいずれかの簡単な簡素化が表示されません。

あなたは

y <- sapply(x, function(u) 
    sum( 
    c(0.5,0.5) 
    %*% 
    mtx.exp(P %*% matrix(c(u,0,0,u),nrow=2), 5) 
) 
) 

ような何かを試みることができるが、あなたのベクトルは本当にめちゃくちゃ長い時間がかかります12^11エントリ、 を持っている場合。 マキシマ、セージ、Yacas：あなたは（2 * 2）の行列、非常に小さな、非常に多数の を持っているので、 あなたが明示的に製品P %*% C といくつかの数式処理システムを使用して、その第五の電力を（計算することができますまた

、 、Mapleなど） と計算された公式を使用してください： これは、ベクトルに対する単純な操作（50行）です。

私はコピーしてRに結果を貼り付け

/* Maxima code */ 
p: matrix([p11,p12], [p21,p22]); 
c: matrix([c1,0],[0,c2]); 
display2d: false; 
factor(p.c . p.c . p.c . p.c . p.c); 

：申し訳ありませんが、多分私は問題を正しく説明していない

c1 <- dnorm(abs(x),0,1); # C is still a diagonal matrix 
c2 <- dnorm(abs(x),1,3); 
p11 <- P[1,1] 
p12 <- P[1,2] 
p21 <- P[2,1] 
p22 <- P[2,2] 
# Result of the Maxima computations: 
# I just add all the elements of the resulting 2*2 matrix, 
# but you may want to do something slightly different with them. 

      c1*(c2^4*p12*p21*p22^3+2*c1*c2^3*p11*p12*p21*p22^2 
           +2*c1*c2^3*p12^2*p21^2*p22 
           +3*c1^2*c2^2*p11^2*p12*p21*p22 
           +3*c1^2*c2^2*p11*p12^2*p21^2 
           +4*c1^3*c2*p11^3*p12*p21+c1^4*p11^5) 
      + 
      c2*p12 
      *(c2^4*p22^4+c1*c2^3*p11*p22^3+3*c1*c2^3*p12*p21*p22^2 
         +c1^2*c2^2*p11^2*p22^2+4*c1^2*c2^2*p11*p12*p21*p22 
         +c1^3*c2*p11^3*p22+c1^2*c2^2*p12^2*p21^2 
         +3*c1^3*c2*p11^2*p12*p21+c1^4*p11^4) 
     + 
     c1*p21 
      *(c2^4*p22^4+c1*c2^3*p11*p22^3+3*c1*c2^3*p12*p21*p22^2 
         +c1^2*c2^2*p11^2*p22^2+4*c1^2*c2^2*p11*p12*p21*p22 
         +c1^3*c2*p11^3*p22+c1^2*c2^2*p12^2*p21^2 
         +3*c1^3*c2*p11^2*p12*p21+c1^4*p11^4) 
     + 
     c2*(c2^4*p22^5+4*c1*c2^3*p12*p21*p22^3 
         +3*c1^2*c2^2*p11*p12*p21*p22^2 
         +3*c1^2*c2^2*p12^2*p21^2*p22 
         +2*c1^3*c2*p11^2*p12*p21*p22 
         +2*c1^3*c2*p11*p12^2*p21^2+c1^4*p11^3*p12*p21)