How to generate samples from a cauchy distribution

For an assignment I have to draw samples of different sizes from a Cauchy distribution in order to show that as the sample size increases, the mean does not tend to any number (therefore the mean does not exist). However, I have been having trouble trying to figure out how to properly create my Cauchy distribution.

Knowing that a cauchy distribution is simply a variable computed using the ratio of two independent standard normal variables, I created two standard normal variables and then divided them to create this new Cauchy variable. I then found the sample mean, and repeated the process for sample sizes of 10, 100, 500, 1000, 5000, and 10,000.

The following code is the same process that I used for each sample size:

a_c1= randn(10,1)
a_c2= randn(10,1)
a_c1ovrc2=a_c1/a_c2
a_c1ovrc2( :, ~any(a_c1ovrc2,1) )= [] %in order to delete all the unnecessary zero columns(I don't understand why they were created in the first place)
a_c1ovrc2_m = mean(a_c1ovrc2)

The above table are the values i obtained after completing each operation separately. I don't feel like this accurately shows that the Cauchy distribution has no mean. Any advice as to what I might be doing wrong would be sincerely appreciated!

Answer 1

[dpb]

a_c1ovrc2( :, ~any(a_c1ovrc2,1) )= []

"in order to delete all the unnecessary zero columns(I don't understand why they were created in the first place)"

They came from the preceding line

a_c1ovrc2=a_c1/a_c2;

HINT: See

doc ./  % follow links to "Arithmetic Operators" at bottom

After perusing these, all should become clear...

[Erol]

I am an idiot and you, Sir, are a genius. Very much appreciated! I definitely need to up my matrix arithmetic skills.

[dpb]

That's probably a little strong on both counts... :)

Answer 2

[John D'Errico]

I recall one simple way to generate Cauchy variates is using the CDF of the Cauchy.

https://en.wikipedia.org/wiki/Cauchy_distribution

The CDF has the form

1/pi*atan(x) + 1 / 2

So if you choose a (uniform) random variate from [0,1], then form the inverse, voila! You have Cauchy deviates.

cvar = @(N) tan((rand(1,N) - 0.5)*pi);

mean(cvar(1000000))
ans =
     -0.18167
mean(cvar(1000000))
ans =
       2.5684
mean(cvar(1000000))
ans =
     0.063165
mean(cvar(1000000))
ans =
      -1.0446
mean(cvar(1000000))
ans =
      0.49277
mean(cvar(1000000))
ans =
     -0.19725
mean(cvar(1000000))
ans =
      0.53434

As you can see, it bounces all over the place, nor does this get better for larger sample size.

mean(cvar(10000000))
ans =
     -0.07473
mean(cvar(10000000))
ans =
      0.78902
mean(cvar(10000000))
ans =
      -1.0978
mean(cvar(10000000))
ans =
      -1.0424

[Erol]

That's a far simpler solution to my problem. Very much appreciated and will absolutely use this method from now on!

[John D'Errico]

Note that a ratio of two normal (zero mean) random variables is apparently also a Cauchy random variable. (I had to check on that.)

https://en.wikipedia.org/wiki/Ratio_distribution

cvar2 = @(N) randn(1,N)./randn(1,N);
>> mean(cvar2(100000))
ans =
       1.4023
>> mean(cvar2(100000))
ans =
      0.34317
>> mean(cvar2(100000))
ans =
     -0.27714
>> mean(cvar2(100000))
ans =
      0.23104
>> mean(cvar2(100000))
ans =
     -0.78052
>> mean(cvar2(100000))
ans =
       3.3391
>> mean(cvar2(100000))
ans =
      -1.3788
>> mean(cvar2(100000))
ans =
      -1.2899
>> mean(cvar2(100000))
ans =
     -0.79651
>> mean(cvar2(100000))
ans =
       1.0549