xrms.cpp

#include <iostream>
#include <random>
#include <cmath>
#include <iomanip>

double exact_rms_normal(int N)
{
    // E[R_N] for N iid N(0,1):
    // sqrt(2/N) * Gamma((N+1)/2) / Gamma(N/2)
	return std::sqrt(2.0 / N) *
			std::tgamma(0.5 * (N + 1)) /
			std::tgamma(0.5 * N);
}

int main()
{
	const int maxN = 50;
	const int n_trials = 1000000;  // adjust for speed/accuracy

    // RNG setup
	std::random_device rd;
	std::mt19937_64 rng(rd());

	std::normal_distribution<double> normal(0.0, 1.0);
	std::student_t_distribution<double> tdist(5.0); // v = 5

    // Scale t_5 to unit variance: Var(t_v) = v / (v - 2)
	const double v = 5.0;
	const double t_var = v / (v - 2.0);      // = 5/3
	const double t_scale = std::sqrt(1.0 / t_var);  // = sqrt(3/5)

    // Accumulators for simulated RMS averages
	double sum_rms_normal[maxN + 1];
	double sum_rms_t[maxN + 1];

	for (int i = 0; i <= maxN; ++i) {
		sum_rms_normal[i] = 0.0;
		sum_rms_t[i] = 0.0;
	}

    // Monte Carlo simulation
	for (int trial = 0; trial < n_trials; ++trial)
	{
		double sumsq_normal = 0.0;
		double sumsq_t = 0.0;

		for (int N = 1; N <= maxN; ++N)
		{
			double z = normal(rng);
			sumsq_normal += z * z;

			double t_raw = tdist(rng);
			double t_std = t_raw * t_scale;  // mean 0, variance 1
			sumsq_t += t_std * t_std;

			double rms_normal = std::sqrt(sumsq_normal / N);
			double rms_t = std::sqrt(sumsq_t / N);

			sum_rms_normal[N] += rms_normal;
			sum_rms_t[N] += rms_t;
		}
	}

	std::cout << std::fixed << std::setprecision(6);
	std::cout << "# n_trials = " << n_trials << "\n";
	std::cout << "N  exact_RMS_normal  sim_RMS_normal  sim_RMS_t5_std\n";

	for (int N = 1; N <= maxN; ++N)
	{
		double exact = exact_rms_normal(N);
		double sim_norm = sum_rms_normal[N] / n_trials;
		double sim_t = sum_rms_t[N] / n_trials;

		std::cout << std::setw(2) << N << "  "
				<< exact << "  "
				<< sim_norm << "  "
				<< sim_t << "\n";
	}

	return 0;
}