Simulate a Beta dist from 2 Gammas

[nickg]

C++ evidently does not include a beta random distribution. However, a Beta(A,B) can be simulated from two gammas: GammaA = Gamma(A,1) and GammaB = Gamma(B,1), by deriving Beta(A,B) as GammaA/(GammaA + GammaB). The C++ code below successfully simulates 1000 random values of each of the two gammas, and creates 1000 presumably correct Betas.

However, when the text files of the three are exported to Matlab for histogram display, the Betas conain some values greater than 1. But running a Matlab script (for i = 1:1000, GammaA(i)/(GammaA(i) + GammaB(i) produces the correct set of random Betas. Why does my C++ code not generate this as I intend?

BetaParameters: [5.0, 10.0, 1000]
 // Distributions.cpp : This file contains the 'main' function. Program execution begins and ends there.
//NG started this "project" Feb 2 2020. For testing simulation of randoms from distributions.
//test beta random sim from two random gammas: X = gamma(A,1), Y = gamma(B,1)
//beta(A,B) = X/(X+Y).
// Runs but the BetaSamples include values > 1. But importing G1Samples and G2Samples into Matlab
//G1Samples./(G1Samples+G2Samples) produces the correct Beta distribution! Need to check why this code
//does not produce BetaSamples identical to G1/(G1+G2).

#include <iostream>
#include <fstream>
#include <random>
#include <string>
using namespace std;

int main()
{
    //   string dir("C:\\Users\\NickG2\\Documents\\Cpp\\VSProjects\\NTU_Examples");
    //   ifstream f_params(dir + "SimParams.txt");
    //   ofstream f_out(dir + "samples.txt");

    ifstream f_params("C:\\Users\\NickG2\\Documents\\Cpp\\VSProjects\\Distributions\\BetaParams.txt");
    ofstream f_out1("C:\\Users\\NickG2\\Documents\\Cpp\\VSProjects\\Distributions\\BetaSamples.txt");
    ofstream f_out2("C:\\Users\\NickG2\\Documents\\Cpp\\VSProjects\\Distributions\\G1Samples.txt");
    ofstream f_out3("C:\\Users\\NickG2\\Documents\\Cpp\\VSProjects\\Distributions\\G2Samples.txt");

    if (!f_params || !f_out1 || !f_out2 || !f_out3) {
        cerr << "Problem opening the files.\n";
        return 1;
    }
    double betaA, betaB;
    int nreps;
    if (!(f_params >> betaA >> betaB >> nreps)) {
        cerr << "Can't read mean and/or stddev and/or nreps.\n";
        return 1;
    }

    gamma_distribution<> dist1(betaA, 1.0);
    gamma_distribution<> dist2(betaB, 1.0);
    random_device rd{};
    mt19937 gen(rd());
    for (int i = 0; i < nreps; ++i) {
        f_out1 << (dist1(gen) / (dist1(gen) + dist2(gen))) << '\n';
        f_out2 << dist1(gen) << '\n';
        f_out3 << dist2(gen) << '\n';
    }
}

[lastchance]

How many distinct random numbers do you use here? (Note that the answer is NOT 2).

        f_out1 << (dist1(gen) / (dist1(gen) + dist2(gen))) << '\n';
        f_out2 << dist1(gen) << '\n';
        f_out3 << dist2(gen) << '\n';

I think that you probably meant this loop to be:

    for (int i = 0; i < nreps; ++i) {
        double X = dist1(gen);
        double Y = dist2(gen);
        f_out1 << X / ( X + Y ) << '\n';
        f_out2 << X << '\n';
        f_out3 << Y << '\n';
    }

[nickg]

Thank you again. I'll try your alternative. I did think my code was only simulating the two gamma randoms.