Home > Back-end >  How to initialize distribution from std::random using conditional logic in C
How to initialize distribution from std::random using conditional logic in C

Time:01-05

I am trying to sample random numbers from different distributions based on conditional logic and I am having trouble finding a nice way to go about it.

I have a struct representing some distribution:

struct Distribution {
    std::string name;
    double args[2];
}

The standard normal distribution would thus be represented by:

Distribution normal = {"normal", {0, 1}}

My goal is to; given an array of Distributions, create a few thousand samples from each distribution. This would not be a problem would it not have been for the fact that the different distributions in std::random are of different type. Due to the different type of the distributions I am having trouble with initializing the distribution before sampling from it.

What I would like to write is something that initializes the distribution outside the sampling loop, e.g.:

struct Distribution {
    std::string name;
    double args[2];
}

int main(void) {
    Distribution distrs[2] = {
        {"uniform", {0, 1}},
        {"normal",  {0, 1}}
    };

    int n_samples = 100;
    double samples[200];

    std::random_device rd;
    std::mt19937 mt(rd());

    for (int ix = 0; ix < 2; ix  ) {
        // Selects distribution to sample from depending on content of distrs[ix]
        some_common_type_for_distr sampler; // This does obv. not compile
        std::string distr_name = distrs[ix].name;
        double args[2] = distrs[ix].args;
        if (distr_name == "uniform") {
            sampler = std::uniform_real_distribution<double>(args[0], args[1]);
        } else if (distr_name == "normal") {
            sampler = std::normal_distribution<double>(args[0], args[1]);
        }

        // Samples from the above initialized distribution
        for (int jx = 0; jx < n_samples; ix  ) { // Sampling loop
            // Shifts by ix*n_samples to place samples from second
            // distribution after samples from first distribution
            samples[jx   ix*n_samples] = sampler(mt);
        }
    }

    // Do something with samples.
}

However, as mentioned above, I run into problems when trying to somehow save the different initialized distributions to the same variable (in above case sampler). This results in me needing to re-initialize the distribution inside the sampling loop, e.g.:

struct Distribution {
    std::string name;
    double args[2];
}

int main(void) {
    Distribution distrs[2] = {
        {"uniform", {0, 1}},
        {"normal",  {0, 1}}
    };

    int n_samples = 100;
    double samples[200];

    std::random_device rd;
    std::mt19937 mt(rd());

    for (int ix = 0; ix < 2; ix  ) {
        // Initializes and samples in one go, resulting in re-initialization
        // on every iteration
        for (int jx = 0; jx < n_samples; ix  ) { // Sampling loop
            std::string distr_name = distrs[ix].name;
            double args[2] = distrs[ix].args;
            if (distr_name == "uniform") {
                std::uniform_real_distribution<double> sampler(args[0], args[1]);
                samples[jx   ix*n_samples] = sampler(mt);
            } else if (distr_name == "normal") {
                std::normal_distribution<double> sampler(args[0], args[1]);
                samples[jx   ix*n_samples] = sampler(mt);
            }
        }    
    }

    // Do something with samples.
}

So my question is basically; is there any way to write the program in such a way that the distributions are initialized outside the sampling loop, as in the first main function above?

CodePudding user response:

you can erase the type, for example, use std::function

std::function<double()> GetSampler(const Distribution& d, std::mt19937& mt){
    if (d.name == "uniform") {
        std::uniform_real_distribution<double> sampler(d.args[0], d.args[1]);
        return [sampler,&mt]()mutable{return sampler(mt);};
    } else if (d.name == "normal") {
        std::normal_distribution<double> sampler(d.args[0], d.args[1]);
        return [sampler,&mt]()mutable{return sampler(mt);};
    }
    // handle other case
    return []{return 4;};
}

use

std::function<double()> samplers[2] = {
    GetSampler(distrs[0],mt),
    GetSampler(distrs[1],mt)
};

// Sampling loop

for (int ix = 0; ix < 2; ix  )
for (int jx = 0; jx < n_samples; jx  ) { 
   auto& sampler = samplers[ix];
   samples[jx   ix*n_samples] = sampler();
}

https://godbolt.org/z/KGsaTTKnr

  • Related