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.

Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s