The SciPy Library (not to be confused with the whole SciPy Project, which includes the other mentioned packages) offers additional functionality on top of NumPy (list is not complete)
Offers nearly 100 continuous and discrete distributions, and statistical tests.
import scipy.stats
x = np.linspace(-5, 5, 100)
standard_normal_pdf = scipy.stats.norm.pdf(x)
plt.plot(x, standard_normal_pdf)
[<matplotlib.lines.Line2D at 0x7f04b93037b8>]
Each distribution has at least two shape parameters:
plt.plot(x, scipy.stats.norm.pdf(x), label='Standard normal')
plt.plot(x, scipy.stats.norm.pdf(x, loc=1), label='Shifted')
plt.plot(x, scipy.stats.norm.pdf(x, scale=2.0), label='Scaled')
plt.plot(x, scipy.stats.norm.pdf(x, loc=-1, scale=0.5), label='Shifted and Scaled')
plt.legend()
plt.show()
Each distributions defines several methods (list not complete):
x = np.linspace(-4, 4, 100)
plt.plot(x, scipy.stats.norm.pdf(x), label='PDF')
plt.legend(loc='upper left')
plt.grid(True)
plt.plot(x, scipy.stats.norm.cdf(x), label='CDF')
plt.legend(loc='upper left')
plt.grid(True)
plt.plot(x, scipy.stats.norm.ppf(x), label='PPF')
plt.legend(loc='upper left')
plt.grid(True)
plt.hist(scipy.stats.norm.rvs(size=1000), bins=20, label='RVS', alpha=0.75)
plt.legend(loc='upper left')
plt.grid(True)
Here's a more advanced example.
We create 1000 samples, each sample contains 20 randomly drawn numbers. Now we perform a statisical test, and check if the 20 randomly drawn numbers are compatible with a gaussian distribution. We do this 1000 times and plot a histogram of the p-values returned by the test.
random_gauss_samples = scipy.stats.norm.rvs(size=(20, 1000))
random_cauchy_samples = scipy.stats.cauchy.rvs(size=(20, 1000))
plt.title('Distribution of P-Values of Normaltest for differently distributed samples')
plt.hist(scipy.stats.normaltest(random_gauss_samples, axis=0).pvalue, bins=10, range=(0,1), label='Gauss', alpha=0.5)
plt.hist(scipy.stats.normaltest(random_cauchy_samples, axis=0).pvalue, bins=10, range=(0,1), label='Cauchy', alpha=0.5)
plt.legend()
plt.show()
Offers:
import scipy.optimize
A transcedental equation like $exp(-x) = -x$, cannot be solved analytically.
One can transform cast an equation into a root finding problem: $$exp(-x) - x = 0$$
scipy.optimize.newton is one possible algorithms which can calculate the root of a function, and therefore solve the equation.
def transcedental_equation(x):
return np.exp(-x) - x
root = scipy.optimize.newton(transcedental_equation, 0.0)
print(root)
0.56714329041
x = np.linspace(-2, 2, 100)
plt.plot(x, np.exp(-x), label='exp(-x)')
plt.plot(x, x, label='x')
plt.plot(x, transcedental_equation(x), label='exp(-x) - x')
plt.axhline(0, color='black', linestyle='--')
plt.ylim((-2, 4))
plt.annotate('Intersection', xy=(root, root), xytext=(1, 3), arrowprops=dict(fc='black', width=0.5, headwidth=5.0))
plt.show()
scipy.special contains the hermite polynomials $H_n$ and the gamma functions $\Gamma(n+1) = n!$.
The wave function of the quantum harmonic oscillator is given by:
$$\varphi(x) = \frac{1}{2^n n!}\ \frac{1}{\pi^{1/4}}\ H_n(x)\ \mathrm{e}^{-x^2}$$import scipy.special
x = np.linspace(-3, 3, 100)
plt.plot(x, x**2, color='black')
for n in range(9):
y = (np.pi**-0.25
* scipy.special.hermite(n)(x) * np.exp(-x**2)
/ (2**n * scipy.special.gamma(n+1))**0.5 )
energy = n + 0.5
plt.fill_between(x, y + energy, energy)
Calculate the probability that a value drawn from a standard normal distribution is in the interval (-1, 1), (-2, 2) and (-3, 3)
Plot the radial component of the wave function for the electron in the hydrogen atom for different quantum numbers $n$ and $l$.
$$R_{nl}(r) = \sqrt{\frac{(n-l-1)!}{2n^4(n+l)!}}\ \mathrm{e}^{-\frac{r}{2n}}\ \left(\frac{r}{n}\right)^l\ L^{2l+1}_{n-l-1}\left(\frac{r}{n}\right)$$$$E_n \sim - \frac{1}{n^2}$$Visualize the Rosenbrock function in a contour-plot and find its minimum
$$ f(x, y) = (1 - x)^2 + 100\: (y - x^2)^2$$