Online Scipy Compiler

SciPy is a collection of numerical routines built on NumPy. Its submodules cover linear algebra (scipy.linalg), optimization (scipy.optimize), integration (scipy.integrate), statistics (scipy.stats), signal processing, interpolation and sparse matrices. Scientists, engineers and students use it to solve equations, fit models to measurements and test whether two groups of data really differ. This page is an online SciPy compiler: the code runs in your browser, so you can try it without installing anything. Run the example first, then paste any snippet below into a new cell to try it.

What the example does

The equations 2x + 3y = 5 and 4x - y = 3 can be written as one matrix equation, A times x, y equals B. Each row of A holds the coefficients of one equation, and B holds the right-hand sides. solve(A, B) from scipy.linalg returns the unknowns as an array, here x = 1.0 and y = 1.0. The result is an array of floats, although every input is an integer. The same call handles more unknowns, as long as A is square: one row per equation and one column per unknown.

Find the minimum of a function

minimize() searches for the inputs that make a function as small as possible, starting from the guess in x0. This function is smallest at x = 3, y = -1:

from scipy.optimize import minimize

def cost(p):
    x, y = p
    return (x - 3) ** 2 + (y + 1) ** 2 + 2

result = minimize(cost, x0=[0, 0])
print(result.success, result.x.round(4), round(result.fun, 4))

result.x holds the best inputs and result.fun the value there. The search stops close to the answer, not exactly on it, which is why the snippet rounds. Check result.success before you trust the result.

Fit a curve to data

curve_fit() finds the parameters of a model function that best match measured data. Here the data comes from an exponential decay with a = 10 and k = 0.8, plus seeded random noise:

import numpy as np
from scipy.optimize import curve_fit

def decay(t, a, k):
    return a * np.exp(-k * t)

rng = np.random.default_rng(seed=0)
t = np.linspace(0, 5, 30)
y = decay(t, 10, 0.8) + rng.normal(0, 0.2, t.size)   # noisy measurements

params, cov = curve_fit(decay, t, y, p0=[5, 1])
errors = np.sqrt(np.diag(cov))
print("a =", params[0].round(3), "+/-", errors[0].round(3))
print("k =", params[1].round(3), "+/-", errors[1].round(3))

p0 is the starting guess. The square roots of the diagonal of cov are the standard errors of the fitted values.

Compute an integral

quad() integrates a function of one variable between two limits, which can be infinite. It returns two numbers: the value and an estimate of its absolute error.

import numpy as np
from scipy.integrate import quad

value, error = quad(np.sin, 0, np.pi)
print(value, error)

# The area under a standard normal curve is 1
area, _ = quad(lambda x: np.exp(-x**2 / 2) / np.sqrt(2 * np.pi), -np.inf, np.inf)
print(round(area, 10))

Run a statistical test

scipy.stats has probability distributions and the common tests. This compares two groups of seeded random measurements with a t-test:

import numpy as np
from scipy import stats

rng = np.random.default_rng(seed=1)
control = rng.normal(loc=50, scale=5, size=40)
treated = rng.normal(loc=53, scale=5, size=40)

result = stats.ttest_ind(control, treated)
print(f"t = {result.statistic:.2f}, p = {result.pvalue:.4f}")

print(stats.norm.cdf(1.96))   # share of a normal distribution below 1.96

For regression models with full summary tables, see statsmodels, which builds on SciPy.

Good to know

• solve() needs a square matrix with exactly one solution. Two equations that describe the same line make the matrix singular, and solve() raises LinAlgError: Matrix is singular. With more equations than unknowns, use scipy.linalg.lstsq() for the least-squares answer.
• curve_fit() starts every parameter at 1 unless you pass p0. For a model such as a sine wave, a poor start can return a wrong fit without any error, so pass a rough guess.
• ttest_ind() assumes both groups have the same variance. Pass equal_var=False for Welch's t-test, which does not.