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

Question

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

我可以找到以下多项式的区域如下:

pn =

   -0.0250    0.0667    0.2500   -0.6000         0

首先使用Simpson规则的整合

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我将使用gammaincie不完整的gamma函数代替gamma)

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

Answer 1

有几件事可能会做.最简单的是避免调用阶乘.相反,人们可以使用那种关系

factorial(n) = gamma(n+1)

由于gamma似乎实际上比对factorial的调用更快,所以你可以节省一点.更好的是,你可以

>> 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个值的时间比计算值多一点.

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