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 Then 使用以下代码使用 gamma 函数评估上述公式中的总和

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 变换从这个分布中获取样本。（为了构建 cdf，我将使用gammainc不完整的 gamma 函数而不是gamma）

我需要从可能具有不同多项式系数的密度中进行采样，并且速度是我关心的问题。我已经可以使用蒙特卡洛方法从这样的密度中采样，但我想看看我是否可以使用密度中的精确采样来加快速度。非常感谢您提前。

4

2 回答 2

7

一个人可以做几件事。最简单的是避免调用阶乘。相反，可以使用以下关系

factorial(n) = gamma(n+1)

由于 gamma 似乎实际上比调用阶乘更快，因此您可以在那里节省一点。更好的是，你可以

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

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

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

更好的是，一个人可以在一次调用 gammaln 时完成所有 4 次调用。例如，想想这是做什么的：

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

请注意，如果您的数字足够大，此调用不会出现溢出问题。而且由于 gammln 是矢量化的，所以一个调用很快。计算 4 个值所花费的时间比计算 1 个值要多一点。

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

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

不可否认，如果你使用 gammaln，你需要在最后调用 exp 来恢复最终结果。但是，您也可以通过一次调用 gamma 来做到这一点。或许是这样的：

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 倍。但它仍然有一组嵌套循环。通过一些努力，您也可以改进它。

例如，与其多次计算所有这些阶乘（或伽马函数），不如执行一次。这应该有效：

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 秒，加速大约为 150 比 1，尽管有三重嵌套循环。

于 2012-11-21T02:19:43.860 回答
4

好吧，甚至无需更改您的代码，您就可以安装来自 Microsoft 的 Tom Minka 的名为lightspeed的出色软件包，它用更快的版本替换了一些内置的 matlab 函数。我知道有一个替代品gammaln()。

您将获得非平凡的速度改进，尽管我不确定多少，而且安装起来很简单。

于 2012-11-21T02:14:10.493 回答