使用 gamma 函数实现嵌套循环的更快方法

Question

我正在尝试评估以下积分：

我可以找到以下多项式的面积如下：

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

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

Answer 1

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

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，尽管有三重嵌套循环。

Answer 2

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

您将获得非同寻常的速度提升，但我不确定提升多少，而且安装起来也很简单。