3

私は次の積分を評価しようとしています：

ここに画像の説明を入力してください

次の多項式の面積は次のように見つけることができます。

pn =

   -0.0250    0.0667    0.2500   -0.6000         0

最初にシンプソンの法則による積分を使用する

fn=@(x) exp(polyval(pn,x));

area=quad(fn,-10,10);
fprintf('area evaluated by Simpsons rule : %f \n',area)

結果はarea evaluated by Simpsons rule : 11.483072 次のようになります。次のコードは、ガンマ関数を使用して上記の式の合計を評価します。

a=pn(1);b=pn(2);c=pn(3);d=pn(4);f=pn(5);
area=0;
result=0;
for n=0:40;
    for m=0:40;
        for p=0:40;
            if(rem(n+p,2)==0)
                result=result+ (b^n * c^m * d^p) / ( factorial(n)*factorial(m)*factorial(p) ) *...
                    gamma( (3*n+2*m+p+1)/4 ) / (-a)^( (3*n+2*m+p+1)/4 );
            end
        end
    end
end

result=result*1/2*exp(f)

そしてこれは11.4831を返します。関数とほぼ同じ結果quad。ここで私の質問は、累積分布関数を作成して、逆CDF変換を使用してこの分布からサンプルを取得できるようにするため、このネストされたループを取り除くことができるかどうかです。（累積分布関数を作成するためにgammainc、代わりに不完全ガンマ関数を使用しますgamma）

多項式係数が異なる可能性のある密度からサンプリングする必要があり、速度が問題になります。私はすでにモンテカルロ法を使用してそのような密度からサンプリングすることができますが、スピードアップするために密度からの正確なサンプリングを使用できるかどうかを確認したいと思います。事前にどうもありがとうございました。

4

2 に答える 2

7

人がするかもしれないいくつかのことがあります。最も簡単なのは、階乗の呼び出しを避けることです。代わりに、次の関係を使用できます

factorial(n) = gamma(n+1)

ガンマは実際には階乗の呼び出しよりも速いように見えるので、そこで少し節約できます。さらに良いことに、あなたはできます

>> timeit(@() factorial(40))
ans =
      4.28681157826087e-05

>> timeit(@() gamma(41))
ans =
      2.06671024634146e-05

>> timeit(@() gammaln(41))
ans =
      2.17632543333333e-05

さらに良いことに、gammalnへの1回の呼び出しで4つの呼び出しすべてを実行できます。たとえば、これが何をするかを考えてみてください。

gammaln([(3*n+2*m+p+1)/4,n+1,m+1,p+1])*[1 -1 -1 -1]'

番号が十分に大きくなった場合でも、この呼び出しはオーバーフローの問題がないことに注意してください。また、gammlnはベクトル化されているため、その1回の呼び出しは高速です。1つを計算するよりも、4つの値を計算する方が少し時間がかかります。

>> timeit(@() gammaln([15 20 40 30]))
ans =
      2.73937416896552e-05

>> timeit(@() gammaln(40))
ans =
      2.46521943333333e-05

確かに、gammalnを使用する場合は、最終結果を復元するために、最後にexpを呼び出す必要があります。ただし、ガンマへの1回の呼び出しでそれを行うこともできます。おそらくこのように：

g = gamma([(3*n+2*m+p+1)/4,n+1,m+1,p+1]);
g = g(1)/(g(2)*g(3)*g(4));

次に、pの内側のループでより創造的になることができます。完全なループではなく、不要な組み合わせを無視するテストと組み合わせて、これを実行してみませんか？

for p=mod(n,2):2:40

そのステートメントは、とにかく使用されたであろうpの値のみを選択するので、ifステートメントを完全に削除できます。

上記のすべてで、ループの速度が約5倍になると思います。ただし、ネストされたループのセットはまだあります。少し努力すれば、それも改善できるかもしれません。

たとえば、これらすべての階乗（またはガンマ関数）を何度も計算するのではなく、1回実行します。これは機能するはずです：

a=pn(1);b=pn(2);c=pn(3);d=pn(4);f=pn(5);
area=0;
result=0;
nlim = 40;
facts = factorial(0:nlim);
gammas = gamma((0:(6*nlim+1))/4);
for n=0:nlim
  for m=0:nlim
    for p=mod(n,2):2:nlim
      result = result + (b.^n * c.^m * d.^p) ...
         .*gammas(3*n+2*m+p+1 + 1) ...
         ./ (facts(n+1).*facts(m+1).*facts(p+1)) ...
         ./ (-a)^( (3*n+2*m+p+1)/4 );
    end
  end
end

result=result*1/2*exp(f)

私のマシンでのテストでは、三重にネストされたループの実行に4.3秒かかることがわかりました。上記の私のバージョンでは同じ結果が得られますが、必要なのは0.028418秒で、ループが3重にネストされているにもかかわらず、約150対1のスピードアップです。

于 2012-11-21T02:19:43.860 に答える
4

コードを変更しなくても、Microsoft の Tom Minka からlightspeedという優れたパッケージをインストールできます。このパッケージは、いくつかの組み込みの matlab 関数をはるかに高速なバージョンに置き換えます。に代わるものがあることは知っていgammaln()ます。

どれくらいかはわかりませんが、自明ではない速度の改善が得られ、インストールも簡単です。

于 2012-11-21T02:14:10.493 に答える