2

I was working with histograms and I have a histogram which has a shape that is well approximated by $a(\cos(x)+1)$. I was wondering if I took a function on the data, say $\cos(x)$ and plotted the resulting data in a histogram, what function would well approximate that histogram? I have written a little code in python to make my question more clear:

import numpy as np
import matplotlib.pyplot as plt

list_of_x = []
list_of_cosx = []


for i in range(100000):
    v = 2np.arccos(2np.random.random()-1)-np.pi
    list_of_x.append(v)
    list_of_cosx.append(np.cos(v))


plt.hist(list_of_x, bins=100)
plt.show()


Well approximated by cos


plt.hist(list_of_cosx, bins=100)
plt.show()


Well approximated by ???


```
lojle
  • 140
Tom
  • 21

1 Answers1

1

When you make a histogram as you did in your python code, what you're actually doing is setting up a probability density function $f_X(x)$ for your variable $x$, then you create a big number $N$ of random numbers in accordance with this probability distribution, and then you count number of occurrences of each possible value of $x$ and you plot it on histogram.

If your familiar with the concept of probability density functions (PDFs), you will notice straight away that this histogram is actually approximating your set PDF, and if you let $N \to \infty$ and the width of the bins (let's call it $\delta$) $\delta \to 0$, you will get the PDF, up to normalization.

Let me put a disclaimer here: I am a complete layman in this field, but I did some research) There is a result in probability theory known as change of variables in PDFs and it says the following. Say you have a random variable $X$ or in other words a PDF $f_X(x)$ defined on reals, and a monotonic function $g: {\mathbb R} \rightarrow {\mathbb R}$ and $g$ is invertible. Then the PDF of a new random variable $Y = g(X)$ is given by: $$f_Y(y) = f_X(g^{-1}(y)) \left\vert \frac d{dy} (g^{-1}(y)) \right\vert$$

Now notice that using this formula we can calculate the PDF for your $x$ values exactly. In your program you used imported random function to generate a random number in interval $[0,1]$ with uniform distribution. We can call this our random $X$ variable with PDF being $$f_X(x) = \left\{\begin{array}{rl} 1, & \text{if } x \in [0,1] \\ 0, & \text{elsewhere } \end{array}\right.$$

And we seek PDF for $Y = g(X)$, where $$g(x) = 2\arccos(2x-1) - \pi$$

We see that $\arccos(x)$ is defined on interval $[-1,1]$ and it is a strictly decreasing invertible function, hence $g(x)$ is too. So we can apply the theorem. First we find that $g$ takes $x \in [0,1]$ (which is just what we want) and the image is interval $[-\pi, \pi]$. So whatever function $f_Y(y)$ we get as PDF of $Y$, we should only look at its part defined on interval $[-\pi,\pi]$. After finding the inverse function $g^{-1}(y)$ and plugging it into the formula provided by the theorem, we get: $$f_Y(y) = \frac14 \sin\left( \frac{y+\pi}2 \right), \quad y \in [-\pi,\pi]$$

Graph of <span class=$f_Y(y)$" />

Now let's compare this with our python generated histogram. To get the width of each bin of the histogram we can say:

$$ \delta = \frac{\text{width of function domain}}{\text{number of bins}} = \frac{2\pi}{100}$$

We will use $f_Y(y)$ to predict the approximate height of a bin centered at $y=0$. The probability of $\vert y \vert \lt \frac{\delta}2$ is $$P(-\frac\delta2 \lt y \lt \frac\delta2) = \int_{-\frac\delta2}^{\frac\delta2} f_Y(y) dy \approx \delta \cdot f_Y(0) = \frac{2\pi}{100} \cdot 0.25$$

This means that if we take $N=100000$ random numbers and plot them, the most likely count of them ending up in this bin is $N$ times the probability we just calculated. If you plug it all in you get that you should expect the peak of your histogram to be somewhere around $1570$. You can run your program and make sure that this is actually the case.

Now, all this partially answers your question. If you take a monotone function of your data you are essentially defining $Z=h(Y)$ and repeating the process.

But if you define $Z=Y^2$, you must be careful to notice that now $y \in [-\pi,\pi]$, but $h(y) = y^2$ is not monotone on this interval. Luckily there is a generalization of the theorem which says that if function $g$ is not monotone on the whole interval i.e. it does not have one unique inverse function, you can essentially divide up that interval into $n$ subintervals on each of which the function is monotone and has an inverse. Then the general formula is:

$$f_Y(y) = \sum_{k=1}^{n} \left| \frac{d}{dy} g^{-1}_{k}(y) \right| \cdot f_X\big(g^{-1}_{k}(y)\big),$$

where $g^{-1}_k(y)$ is an inverse of $g(x)$ on the $k$-th such subinterval.

In our example we divide $[-\pi, \pi]$ into $[-\pi,0]$ and $(0,\pi]$. $$h(y) = y^2 \quad \implies \quad h^{-1}_{1,2}(z) = \pm \sqrt{z},$$ where $+\sqrt{z}$ is inverse on $(0,\pi]$ and $-\sqrt{z}$ is inverse on $[-\pi,0]$. After doing the dirty work, you end up with a result which you can plot and check again with python using the same technique described above.

$$f_Z(z) = \frac1{8\sqrt{z}} \left[ \sin\left(\frac{\pi+\sqrt{z}}{2}\right) + \sin\left(\frac{\pi-\sqrt{z}}{2}\right) \right], \quad z \in [0,\pi^2]$$

lojle
  • 140