Sampling from Cauchy (and Gaussian) distribution

How we construct  a random draw  from a Cauchy distribution, if we start with a draw from uniform distribution? That’s what asks the Exercise 7.3 of the book The Concepts and Practice of Mathematical Finance by Mark Joshi. We post here some C++ code as an answer. This is easily done by explicitly computing the inverse cumulative distribution function F^{-1}(x) of the Cauchy distribution and plugging then random numbers generated with a uniform distribution (this is done with the function GetOneRandomByICD).

f(x) = \frac{1}{\pi}\frac{1}{1+x^2},\; F(x) = \frac{1}{\pi}\arctan (x) + \frac{1}{2},
F^{-1}(x) = \tan [\pi (x-1/2)]

More tricky is the case in which we don’t have the F^{-1}(x). In that case we can obtain our sample by using two random uniforms – GetOneCauchyByUniforms –  or using two random numbers with gaussian distribution – GetOneCauchyByGaussians. In the latter case of course we need also a routine  to compute gaussian numbers. We present here the routine GetOneGaussianByRatio, which is similar to GetOneCauchyByUniforms (both are based on ratio of uniforms – see Probability and Random Processes by Grimmet and Stirzaker). The mathematical proofs are in the probability section of my notes.

#ifndef RANDOM_H
#define RANDOM_H

double CauchyICD(double uniform);
double GetOneRandomByICD(double (*func_ptr)(double));

double GetOneCauchyByUniforms();
double GetOneCauchyByGaussian();

double GetOneGaussianByRatio();

#include "Random.h"
#include <cmath>
#include <cstdlib>

double CauchyICD(double x){
 return std::tan(M_PI*(x-0.5));
double GetOneRandomByICD(double (*func_ptr)(double)){
	return (*func_ptr)(rand()/static_cast<double>(RAND_MAX));

double GetOneGaussianByRatio(){

	double x;
	double y;

	double ratio;

	do {
		x = 2.*rand()/static_cast(RAND_MAX)-1.;
		y = 2.*rand()/static_cast(RAND_MAX)-1.;
		ratio = x/y;
	while((y > std::exp(-ratio*ratio*0.25))     
           || (y < -std::exp(-ratio*ratio*0.25)) );

	return ratio;

double GetOneCauchyByUniforms(){

	double x;
	double y;

	double sizeSquared;

	do {
		x = 2.*rand()/static_cast<double>(RAND_MAX)-1.;
		y = 2.*rand()/static_cast<double>(RAND_MAX)-1.;
		sizeSquared = x*x + y*y;
	while(sizeSquared >= 4.0);

	return y/x;

double GetOneCauchyByGaussian(){
	double x = GetOneGaussianByRatio();
	double y = GetOneGaussianByRatio();

	return x/y;

A simple test program which write to files samples of Cauchy distributed numbers:

#include "Random.h"
#include <iostream>
#include <fstream>
#include <cmath>
#include <cstdlib>

int main(){

	long sampleLength = 10000;

	std::ofstream myfile1; ("inverse_cauchy.txt");
	std::ofstream myfile2; ("ratio_cauchy.txt");	
        std::ofstream myfile3; ("gauss_cauchy.txt");	

    for(long n=0; n< sampleLength; n++){
		 myfile1 << GetOneRandomByICD(&CauchyICD) << "\n";
		 myfile2 << GetOneCauchyByRatio() << "\n";
		 myfile3 << GetOneCauchyByGaussian() << "\n";


	return 0;

Now just compile and try it.