How to generate a Cauchy random variable

In general, the fundamental result you need is that for any random variable $X$ with distribution function $F$, the random variable $Y=F(X)$ has a uniform distribution on $[0,1]$. Consequently, if you can invert $F$, then you can use a uniform density to simulate back your random variable $X$, since $X = F^{-1}(Y)$.'

In your case, the cdf of a standard Cauchy is $$F=\dfrac{1}{\pi}\arctan(x)+\dfrac{1}{2}$$ and therefore if you let $$y= \dfrac{1}{\pi}\arctan(x)+\dfrac{1}{2}$$ you immediately get $$x = \tan{\left(\pi \left(y - \dfrac{1}{2}\right)\right)}$$

Hence, to generate a standardized Cauchy, use the rand function in Matlab to generate a uniform $[0,1]$ variate subtract 1/2 from it, multiply the result by $\pi$, and apply the tangent function. Repeat a bunch of times to get your sample.

Another interesting way to simulate a Cauchy variable is based upon the observation that the ratio of two Normal distributions is Cauchy distributed. Hence you can generate two standard independent Normal arrays and do their ratio term by term. The resulting array will be Cauchy distributed.

The expected value of a random variable with density $f$ is equal to $$ \int_{-\infty}^\infty xf(x)\,dx $$ provided $$ \int_{-\infty}^\infty |x|f(x)\,dx<\infty. $$ If that last condition fails to hold, then some strange things can happen. One of those is that the conclusion of the strong law of large numbers will fail. With the Cauchy distribution, one has $$ \int_{-\infty}^\infty |x|f(x)\,dx=\infty. $$ Consequently things like the following can happen: \begin{align} & \int_{-a}^a x f(x)\,dx \to 0\text{ as } a\to\infty \\[8pt] \text{ but } & \int_{-a}^{2a} xf(x)\,dx\to\text{a nonzero number as }a\to\infty. \end{align} Work out the integrals and you'll see that. That's the "heavy-tailed" nature of the thing.

Now suppose $X$ is uniformly distributed between $-\pi/2$ and $\pi/2$. Then \begin{align} f_{\tan X}(w) & =\frac{d}{dw} F_{\tan X}(w) = \frac{d}{dw}\Pr(\tan X\le w) \\[8pt] & = \frac{d}{dw} \Pr(X\le\arctan w) \\[8pt] & = \frac{d}{dw}\int_{-\pi/2}^{\arctan w} \frac{dx}\pi \\[8pt] & = \frac1\pi\cdot\frac{d}{dz} \left( \arctan w +\frac\pi2 \right) \\[8pt] & = \frac{1}{\pi(1+w^2)}. \end{align}

So there you have a way to simulate a Cauchy-distributed random variable: First simulate a random variable uniformly distributed between $\pm\pi/2$. Then take its tangent (in radians).

Cauchy random variables can be obtained by sampling the inverse CDF of the distribution. This technique is referred to as inverse transform sampling and is very useful for generating random variates from many distributions. Starting with the cumulative distribution function, where $\mu$ is the location parameter and $c$ is the scale parameter

$$F(x;\mu,c) = \frac{1}{\pi}\text{arctan}\left(\frac{x - \mu}{c}\right) - \frac{1}{2}$$

we solve for $x$ as a function for $F$, which in this case is straightforward:

$$x = c\cdot\text{tan}(\pi (F - ½))+\mu$$

$F$, being the CDF, varies from 0 to 1. Thus, if we can replace $F$ with values randomly sampled from the uniform distribution on $(0, 1)$, we obtain corresponding values sampled from the Cauchy distribution. If $\mu = 0$ and $c = 1$ one obtains the "standard" Cauchy distribution.

In Matlab:

function r = cauchyrnd(mu,c,varargin)
%R = CAUCHYRND(mu,c,M,N)
r = c.*tan(pi*(rand(varargin{:})-0.5))+mu;

The first two arguments are the location and scale parameters of the Cauchy distribution and any subsequent arguments just need to match the input format of the rand function.