使用 gsl_histogram_pdf_sample 从临时分布中采样

Question

我在理解如何在 C++ 中使用库 GSL 中的函数 gsl_histogram_pdf_sample 时遇到了一些麻烦。文档在这里，

http://www.gnu.org/software/gsl/manual/html_node/The-histogram-probability-distribution-struct.html#The-histogram-probability-distribution-struct

我还不是真正的专家，所以，我想知道是否有人可以告诉我这段代码有什么问题，

#include <iostream>
#include <gsl/gsl_histogram.h>
#include <gsl/gsl_rng.h>

using namespace std;

int main()
{
    // I am going to use 5 bins
    size_t  Bins = 5;
    // These are the ranges (must be Bins + 1)
    double range[6] = { 0.0, 1.0, 2.0, 3.0, 4.0, 5.0 };
    // Array with probabilities
    double w[5]  = {0.05, 0.1, 0.3, 0.4, 1};

    // Create the histogram pdf

    gsl_histogram_pdf MyHistPdf;

    MyHistPdf.n     = Bins;
    MyHistPdf.range = range;
    MyHistPdf.sum   = w;

    const gsl_rng_type * T;
    gsl_rng * r;
    T = gsl_rng_default;
    r = gsl_rng_alloc (T);
    double u = gsl_rng_uniform(r);
    cout << u << endl;
    double a = gsl_histogram_pdf_sample(&MyHistPdf, u);

    return 0;
}

程序编译没有错误，但运行时总是报如下错误，

gsl: /usr/src/gsl-1.16-1/src/gsl-1.16/histogram/pdf.c:46: ERROR: cannot find r in cumulative pdf

我不知道这是什么意思。

Answer 1

在回答您的问题之前，您需要了解有关 GSL 的基本规则：“除非您详细了解 GSL 源代码，否则切勿直接操作 gsl 结构变量”。如果你这样做，你的程序将不可避免地失败，或者你的结果将毫无意义！为什么？因为您不详细了解有关这些变量如何存储和操作的内部 gsl 约定！

仔细阅读文档，您会发现有 API 函数可以为您分配和操作所有结构变量。如果您是 C++ 编码员，一个好的“经验法则”是想象 GSL 结构中的所有变量都是私有的（不是这种情况的唯一原因是因为 C 没有使成员变量私有的能力！）。

在文档中很清楚需要分配 gsl_histogram 结构，需要插入一些数据，并且您使用的 gsl_histogram_pdf 需要使用特定的 API 函数进行初始化，然后才能调用 gsl_histogram_pdf_sample。我看到你试图通过破解“sum”指针来避免插入一些数据，而 GSL 失败了，正如我所说，如果你不遵循基本规则，就会发生这种情况！

这是我写的一个例子来指导你

#include <iostream>
#include <cassert>
#include <gsl/gsl_histogram.h>
#include <gsl/gsl_rng.h>

using namespace std;

int main()
{
   size_t  Bins = 5;
   double range[6] = { 0.0, 1.0, 2.0, 3.0, 4.0, 5.0 };

    // histogram 
    gsl_histogram* hist = gsl_histogram_alloc (Bins);
    assert( hist != NULL );
    gsl_histogram_set_ranges (hist, range, 6);

    // include some elements
    const gsl_rng_type * T = gsl_rng_default;
    gsl_rng *r = gsl_rng_alloc (T);

    for(int j = 0; j < 10000; ++j ) {
      const double u = gsl_rng_uniform(r) * 5.0; 
      gsl_histogram_increment (hist, u);            
    }

    // Create the histogram pdf
    gsl_histogram_pdf*  MyHistPdf = gsl_histogram_pdf_alloc (Bins);
    assert( MyHistPdf != NULL );
    int status = gsl_histogram_pdf_init (MyHistPdf, hist);
    assert( status != GSL_EDOM );

    std::cout << "result you want";
    std::cout << gsl_histogram_pdf_sample(MyHistPdf,gsl_rng_uniform(r)) << std::endl;           
    return 0;
  }