# 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.

//Random.h
#ifndef RANDOM_H
#define RANDOM_H

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

double GetOneCauchyByUniforms();
double GetOneCauchyByGaussian();

double GetOneGaussianByRatio();

#endif
//Random.cpp
#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:

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

int main(){

long sampleLength = 10000;

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

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

myfile1.close();
myfile2.close();
myfile3.close();

return 0;
}

Now just compile and try it.