草庐IT

c++ - 概率密度的正交例程

coder 2024-02-05 原文

我想对 (-\infty, a] 中的概率密度函数进行积分，因为 cdf 在封闭形式下不可用。但我不确定如何在 C++ 中执行此操作。

这个任务在 Mathematica 中非常简单；我需要做的就是定义函数，

f[x_, lambda_, alpha_, beta_, mu_] := 
   Module[{gamma}, 
     gamma = Sqrt[alpha^2 - beta^2]; 
     (gamma^(2*lambda)/((2*alpha)^(lambda - 1/2)*Sqrt[Pi]*Gamma[lambda]))*
      Abs[x - mu]^(lambda - 1/2)*
      BesselK[lambda - 1/2, alpha Abs[x - mu]] E^(beta (x - mu))
   ];

然后调用NIntegrate例程对其进行数值积分。 

F[x_, lambda_, alpha_, beta_, mu_] := 
    NIntegrate[f[t, lambda, alpha, beta, mu], {t, -\[Infinity], x}]

现在我想用 C++ 实现同样的事情。我使用 gsl 数值库中的例程 gsl_integration_qagil。它旨在集成半无限区间 (-\infty, a] 上的函数，这正是我想要的。但不幸的是我无法让它工作。

这是C++中的密度函数，

density(double x)
{
using namespace boost::math;

if(x == _mu)
    return std::numeric_limits<double>::infinity();

    return pow(_gamma, 2*_lambda)/(pow(2*_alpha, _lambda-0.5)*sqrt(_pi)*tgamma(_lambda))* pow(abs(x-_mu), _lambda - 0.5) * cyl_bessel_k(_lambda-0.5, _alpha*abs(x - _mu)) * exp(_beta*(x - _mu));

}

然后我尝试通过调用 gsl 例程进行集成以获取 cdf。 

cdf(double x)
{
gsl_integration_workspace * w = gsl_integration_workspace_alloc (1000);

    double result, error;      
    gsl_function F;
    F.function = &density;

    double epsabs = 0;
    double epsrel = 1e-12;

    gsl_integration_qagil (&F, x, epsabs, epsrel, 1000, w, &result, &error);

    printf("result          = % .18f\n", result);
    printf ("estimated error = % .18f\n", error);
    printf ("intervals =  %d\n", w->size);

    gsl_integration_workspace_free (w);

    return result;

}

但是 gsl_integration_qagil 返回错误，迭代次数不足。 

 double mu = 0.0f;
 double lambda = 3.0f;
 double alpha = 265.0f;
 double beta = -5.0f;

 cout << cdf(0.01) << endl;

如果我增加工作区的大小，则贝塞尔函数将不会计算。

我想知道是否有人可以让我深入了解我的问题。使用 x = 0.01 调用上面相应的 Mathematica 函数 F 返回 0.904384。

会不会是密度集中在一个非常小的区间附近(即在[-0.05, 0.05] 之外密度几乎是0，给出了一个图以下)。如果是这样，可以做些什么。谢谢阅读。

最佳答案

回复:积分到 +/- 无穷大:

我会使用 Mathematica 找到 |x - μ| 的经验界限>> K，其中 K 表示均值周围的“宽度”，K 是 alpha、beta 和 lambda 的函数——例如 F 小于且约等于 a(x-μ)-2 或 ae-b(x-μ)2 或其他。这些函数具有无穷大的已知积分，您可以根据经验对其进行评估。然后您可以对 K 进行数值积分，并使用有界近似从 K 到无穷大。

计算出 K 可能有点棘手；我对贝塞尔函数不是很熟悉，所以我帮不了你太多。

一般来说，我发现对于不明显的数值计算，最好的方法是在进行数值评估之前尽可能多地进行分析数学运算。 (有点像自动对焦相机 - 将其靠近您想要的位置，然后让相机完成剩下的工作。)

关于c++ - 概率密度的正交例程，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/10776805/

amp43codelambdadoublec++mathwolfram-mathematica

有关c++ - 概率密度的正交例程的更多相关文章

  1. ruby-on-rails - 如何优雅地重启 thin + nginx？ - 2

    我的瘦服务器配置了nginx，我的ROR应用程序正在它们上运行。在我发布代码更新时运行thinrestart会给我的应用程序带来一些停机时间。我试图弄清楚如何优雅地重启正在运行的Thin实例，但找不到好的解决方案。有没有人能做到这一点？ 最佳答案 #Restartjustthethinserverdescribedbythatconfigsudothin-C/etc/thin/mysite.ymlrestartNginx将继续运行并代理请求。如果您将Nginx设置为使用多个上游服务器，例如server{listen80;server

  2. ruby - 使用 `+=` 和 `send` 方法 - 2

    如何将send与+=一起使用？a=20;a.send"+=",10undefinedmethod`+='for20:Fixnuma=20;a+=10=>30 最佳答案 恐怕你不能。+=不是方法，而是语法糖。参见http://www.ruby-doc.org/docs/ProgrammingRuby/html/tut_expressions.html它说Incommonwithmanyotherlanguages,Rubyhasasyntacticshortcut:a=a+2maybewrittenasa+=2.你能做的最好的事情是:

  3. ruby - 如何计算 Liquid 中的变量 +1 - 2

    我对如何计算通过{%assignvar=0%}赋值的变量加一完全感到困惑。这应该是最简单的任务。到目前为止，这是我尝试过的:{%assignamount=0%}{%forvariantinproduct.variants%}{%assignamount=amount+1%}{%endfor%}Amount:{{amount}}结果总是0。也许我忽略了一些明显的东西。也许有更好的方法。我想要存档的只是获取运行的迭代次数。 最佳答案 因为{{incrementamount}}将输出您的变量值并且不会影响{%assign%}定义的变量，我

  4. arrays - Ruby 数组 += vs 推送 - 2

    我有一个数组数组，想将元素附加到子数组。+=做我想做的，但我想了解为什么push不做。我期望的行为(并与+=一起工作):b=Array.new(3,[])b[0]+=["apple"]b[1]+=["orange"]b[2]+=["frog"]b=>[["苹果"],["橙子"],["Frog"]]通过推送，我将推送的元素附加到每个子数组(为什么？):a=Array.new(3,[])a[0].push("apple")a[1].push("orange")a[2].push("frog")a=>[[“苹果”、“橙子”、“Frog”]、[“苹果”、“橙子”、“Frog”]、[“苹果”、“

  5. += 的 Ruby 方法 - 2

    有没有办法让Ruby能够做这样的事情？classPlane@moved=0@x=0defx+=(v)#thisiserror@x+=v@moved+=1enddefto_s"moved#{@moved}times,currentxis#{@x}"endendplane=Plane.newplane.x+=5plane.x+=10putsplane.to_s#moved2times,currentxis15 最佳答案 您不能在Ruby中覆盖复合赋值运算符。任务在内部处理。您应该覆盖+，而不是+=。plane.a+=b与plane.a=

  6. ruby - Sinatra + Heroku + Datamapper 使用 dm-sqlite-adapter 部署问题 - 2

    出于某种原因，heroku尝试要求dm-sqlite-adapter，即使它应该在这里使用Postgres。请注意，这发生在我打开任何URL时-而不是在gitpush本身期间。我构建了一个默认的Facebook应用程序。gem文件:source:gemcuttergem"foreman"gem"sinatra"gem"mogli"gem"json"gem"httparty"gem"thin"gem"data_mapper"gem"heroku"group:productiondogem"pg"gem"dm-postgres-adapter"endgroup:development,:t

  7. ruby - Ruby 中字符串运算符 + 和 << 的区别 - 2

    我是Ruby和这个网站的新手。下面两个函数是不同的，一个在函数外修改变量，一个不修改。defm1(x)x我想确保我理解正确-当调用m1时，对str的引用被复制并传递给将其视为x的函数。运算符当调用m2时，对str的引用被复制并传递给将其视为x的函数。运算符+创建一个新字符串，赋值x=x+"4"只是将x重定向到新字符串，而原始str变量保持不变。对吧？谢谢 最佳答案 String#+::str+other_str→new_strConcatenation—ReturnsanewStringcontainingother_strconc

  8. ruby - rails 3.2.2(或 3.2.1)+ Postgresql 9.1.3 + Ubuntu 11.10 连接错误 - 2

    我正在使用PostgreSQL9.1.3(x86_64-pc-linux-gnu上的PostgreSQL9.1.3，由gcc-4.6.real(Ubuntu/Linaro4.6.1-9ubuntu3)4.6.1，64位编译)和在ubuntu11.10上运行3.2.2或3.2.1。现在，我可以使用以下命令连接PostgreSQLsupostgres输入密码我可以看到postgres=#我将以下详细信息放在我的config/database.yml中并执行“railsdb”，它工作正常。开发:adapter:postgresqlencoding:utf8reconnect:falsedat

  9. ruby - 在 Ruby + Chef 中检查现有目录失败 - 2

    这是我在ChefRecipe中的一blockRuby:#ifdatadirdoesn'texist,moveoverthedefaultoneif!File.exist?("/vol/postgres/data")execute"mv/var/lib/postgresql/9.1/main/vol/postgres/data"end结果是:Executingmv/var/lib/postgresql/9.1/main/vol/postgres/datamv:inter-devicemovefailed:`/var/lib/postgresql/9.1/main'to`/vol/post

  10. ruby-on-rails - 使用 Pow 作为服务器在 RubyMine 中调试 - Ruby 2.1.1 + Rails 4 - 2

    我已经开始使用RubyMine6。我正在处理Rails4、Ruby2.1.1项目。我无法找到如何使用Pow作为服务器调试到RubyMine。你能给我指明正确的方向吗？ 最佳答案 我能够使用远程调试从RubyMine进行调试。我正在使用RubyMine6、Rails3、Ruby2.1.1。首先创建一个.powenv文件并添加:exportRUBY_DEBUG_PORT=1234exportPOW_WORKERS=1将以下gem添加到您的Gemfile:gem'ruby-debug-ide'gem'debase'创建一个新的初始化器st

随机推荐