Random variables#
A random variable X is a variable whose value depends on the outcomes of a random phenomenon. \(X(\omega)\) is a function from possible outcomes of a sample space \(\omega \in \Omega\).
For a coin toss \(\Omega={H,T}\) \(X(H)=+1\) and \(X(T)=-1\). Every time the experiment is done, X returns either +1 or -1. We could also make functions of random variables, e.g., every time X=+1, we ear 25 cents, etc.
Random variables are classified into two main types: discrete and continuous.
Discrete Random Variable: It assumes a number of distinct values. Discrete random variables are used to model scenarios where outcomes can be counted, such as the number of particles emitted by a radioactive source in a given time interval or the number of photons hitting a detector in a certain period.
Continuous Random Variable: It can take any value within a continuous range. These variables describe quantities that can vary smoothly, such as the position of a particle in space, the velocity of a molecule in a gas, or the energy levels of an atom.
Probability distribution of random variable#
For any random variable X, we are interested in finding probability distribution over possible values \(x\) that \(X\) takes: \(p_X(x)\).
Note the very subtle but crucial difference between \(x\) (1-6 for a die) and a random variable \(X\), which generates numbers \(x\) according to distributions \(p(x)\)
Expectation#
\(E[x]\) is the theoretical value of the mean as opposed to the sample mean which we compute in simulations. For instance, consider the difference between the average height of people calculated from a small sample of cities versus that from the entire world population. Will see that with enough sampling the sample mean approaches its expectation. expectation can apply to 1 (normalization) to variable (mean) or any function of variable, e.g mean subtracted square (variance).
Variance (fluctuation)#
Variance is an expectation of squared deviation of x from expectation, which characterizes the spread in values of x.
Sum of two random variables#
Consider sum of two random variables, suppose these are sum of numbers from two dies or sum of two coin flips. Obviously sums of random variables are random variables themselves.
Expectation is always a linear operator as you can verify by plugging in definition and using linear property of definite integrals.
We see that Variance, in general, is not a linear operator.
Denoting mean subtracting variables \(Y_i=X_i-E[X_i]\) we can show this by expanding square in the expectation expression.
We will see in the next section that for independent random variables the last cross terms is zero. This property has huge significance for statistical mechanics, statistics and sciences in general
Random numbers in python#
The numpy.random has the fastest random number generators based on low-level code written in C.
The Scipy.stats has an extensive library of statistical distributions and tools for statistical analysis.
First, we take a look at the most widely used random numbers of numpy, also called standard random numbers. These are rand (uniform random number on interval 0,1) and randn (stnadard average random number with 0 mean and 1 variance).
When running code that uses random numbers results will always differ for every run. If you want code to reproduce the same result, you can fix the seed to get reproducible results:
np.random.seed(8376743)To convert random variables to probability distributions we need to generate large enough sample then perform histogramming via
np.histor directly histogram and visualize by one shot viaplt.hist()
import numpy as np
X = np.random.rand(30)
print(X)
plt.plot(X)
counts, edges = np.histogram(X, range=(0,1), bins=20)
counts, edges
plt.hist(X, density=True)
Probability distributions#
def rnplot(r):
'''Convenience function for making quick two-panel plot showing
a line plot for the sequence of random numbers (RN)
a histogram plot of the probability density of random numbers
'''
fig, ax = plt.subplots(ncols=2)
ax[0].plot(r, color='blue', label='trajectory')
ax[1].hist(r, density=True, color='red', label = 'histogram')
ax[0].set_xlabel('Samples of RN')
ax[0].set_ylabel('Values of RN')
ax[1].set_xlabel('Values of RN')
ax[1].set_ylabel('Probability Density')
fig.legend();
fig.tight_layout()
Uniform#
Probability distribution
\(E[x] = \frac{1}{2}(a+b)\)
\(V[x] = \frac{1}{12}(a-b)^2\)
Random Variable
\(U(a, b)\) modeled by
np.random.uniform(a,b, size=(N, M))\(U(0, 1)\) modeled by
np.random.rand(N, M, P, ...)
np.random.rand(5)
r = np.random.rand(200)
rnplot(r)
Binomial#
Probability Mass Function#
\(E[n] = Np\)
\(V[n] = 4Np(1-p)\)
Random Variable
\(B(n, p)\) modeled by
np.random.binomial(n, p, size)
np.random.binomial(n=10, p=0.6, size=20)
r = np.random.binomial(n=10, p=0.6, size=2000)
rnplot(r)
Gaussian#
\(E[x] = \mu\)
\(V[x] = \sigma^2\)
Random Variable
\(N(a, b)\) modeled by
np.random.normal(loc,scale, size=(N, M))\(N(0, 1)\) modeled by
np.random.randn(N, M, P, ...)
# For a standard normal with sigma=1, mu=0
r = np.random.randn(200)
# For a more general Gaussian
rgauss = np.random.normal(loc=2., scale=5, size=200)
rnplot(r)
Transforming random variables#
Changing from random variable X to a new variable \(Y=f(X)\) leads to a relation between two probability functions where the multiplying factor is jacobian. E.g think of a factor accounting for changing the spacing between points (or volumes) \(p(x)dx=p(y)dy\)
Let us consider addition \(Y=X+a\) where we have \(p(y)=p(x+a) \cdot 1\) and multiplication \(Y=aX\) where we have \(p(y) = p(x)\cdot \frac{1}{a}\)
These transformations mean we can shift average by a constant factor \(E[X+a]=E[X]+a\). And when the original variable is multiplied by a constant factor \(a\) the variance become \(a^2\) times more: \(V[ax] = a^2V[x]\).
As an example we can use standard uniform and standard normal to carry out the following transformations:
Exact vs sampled probability distributions#
scipy.stats contains a large number of probability distributions. Explores examples there
from scipy.stats import binom, norm, poisson
xmin, xmax, step = -5, 5, 1000
dx = (xmax-xmin)/step
x = np.linspace(xmin, xmax, step)
px = norm.pdf(x)
print('normalization', np.sum(px*dx))
r = np.random.randn(100)
plt.hist(r, density=True, label=f'Sampled mean={r.mean():.2f}, var={r.var():.2f}');
plt.plot(x, px,'k', label='Exact mean=0, var=1')
plt.legend(loc=2)
plt.ylabel('$p(x)$')
plt.xlabel('$x$')
Problems#
Problem 1 Binomial as generator of Gaussian and Poisson distributions#
Show that in large number limit binomial distribution tends to gaussian. Show is by expanding binomial distirbution \(logp(n)\) in power series showing that terms beyond quadratic can be ignored.
In the limit \(N\rightarrow \infty\) but for very small values of \(p \rightarrow 0\) such that \(\lambda =pN=const\) there is another distribution that better approximates Binomial distribution: \(p(x)=\frac{\lambda^k}{k!}e^{-\lambda} \) It is known as Poisson distribution.
Poisson distribution is an excellent approximation for probabilities of rare events. Such as, infrequently firing neurons in the brain, radioactive decay events of Plutonium or rains in the desert.
Derive Poisson distribution by taking the limit of \(p\rightarrow 0\) in binomial distribution.Using numpy and matplotlib plot binomial probability distribution against Gaussian and Poisson distributions for different values of N=(10,100,1000,10000).
For a value N=10000 do four plots with the following values p=0.0001, 0.001, 0.01, 0.1. You can use subplot functionality to make a pretty 4 column plot. (See plotting module)
fig, ax = plt.subplots(nrows=1, ncols=4)
ax[0].plot()
ax[1].plot()
ax[2].plot()
ax[3].plot()