Symbolic Regression in Python in 3 Steps, with a Real Example

Updated November 8, 2026. Tested with turingbot 4.0 (library and program), Python 3.13, pandas 3.1, scikit-learn 1.9 and sympy 1.14.

Symbolic regression finds a mathematical formula that fits your data, instead of a black-box model. This tutorial runs it from Python with TuringBot's library, on a dataset where we know the right answer, so you can see whether the search finds it. The search takes 30 seconds and runs in the free version, and the result is a formula you can evaluate with NumPy, export to C, or check against theory.

Step 1: Install

The search engine is the TuringBot program, and the Python library runs it in the background. Install both:

1. Download TuringBot for Windows, macOS or Linux and install it.
2. Install the library from PyPI:

pip install turingbot pandas scikit-learn sympy

The library finds the TuringBot installation by itself. If you installed it somewhere unusual, pass its location with start_process(path=...). pandas, scikit-learn and sympy are only for the examples below.

Step 2: Run the search

The data: 200 throws of a ball, with the launch speed, the launch angle, and how far it went. We generate it from the physics formula, distance = speed² · sin(2 · angle) / 9.81, and add 1% of random noise, like a measurement error. The search doesn't know the formula; the question is whether it can find it.

import time

import numpy as np
import pandas as pd
import turingbot as tb
from sklearn.model_selection import train_test_split

# 200 throws: launch speed (m/s), launch angle (radians), distance (m), with 1% noise
rng = np.random.default_rng(0)
df = pd.DataFrame({
    'speed': rng.uniform(5, 30, 200),
    'angle': rng.uniform(0.1, 1.4, 200)
})
df['distance'] = df.speed**2 * np.sin(2 * df.angle) / 9.81 * rng.normal(1, 0.01, 200)
train, test = train_test_split(df, test_size=0.25, random_state=0)

sim = tb.simulation()
sim.start_process(input_file=train)  # The last column is the target
time.sleep(30)  # Let it search for 30 seconds
sim.refresh_functions()
sim.terminate_process()

for complexity, error, formula in sim.functions:
    print(complexity, error, formula)

start_process takes a DataFrame, a NumPy array or a file, and starts the search in the background, so your script or notebook stays free while it runs. refresh_functions() loads the formulas found so far into sim.functions: the best formula at each complexity, where a larger formula appears only if it beats every smaller one. Here are selected rows from our run (constants rounded; the error is the RMS error in meters, on the training rows; complexities 4, 7 and 10, small steps such as lgamma(speed), are left out):

The story is in the error column. The search first finds that distance grows with speed squared (complexity 5), then that it rises and falls with the angle (complexity 11, a parabola in the angle), and at complexity 13 it finds the actual law: 0.10198 * speed² * sin(1.9992 * angle). That is speed² · sin(2 · angle) / 9.81 (1/9.81 = 0.10194), with the constants slightly off because of the noise.

Its error, 0.404 m, is the size of the noise we added (0.405 m on these rows). So no formula can honestly do better. The search went on to complexity 59 in 30 seconds, and those larger formulas lower the training error by about 9 cm (to 0.311 m) by fitting the noise. On the test rows the search never saw, the complexity 59 formula does worse: 0.435 m, against 0.366 m for complexity 13. That is how you choose a formula in general: go down the list to where the error stops falling fast, and take the formula there.

Step 3: Use the formula

The formulas are plain text, so sympy can parse them, turn them into NumPy functions, and score them on the test rows the search never saw:

import sympy
from sklearn.metrics import r2_score

formulas = {complexity: formula for complexity, error, formula in sim.functions}
expr = sympy.sympify(formulas[13])
f = sympy.lambdify(['speed', 'angle'], expr, 'numpy')

print(r2_score(test.distance, f(test.speed, test.angle)))  # 0.9998

sympy reads most formulas as they are. A few of TuringBot's functions need a translation first; complexity 16 above uses round, for example:

names = {
    'round': lambda x: sympy.floor(x + sympy.Rational(1, 2)),
    'ceil': sympy.ceiling,
    'lgamma': sympy.loggamma,
    'tgamma': sympy.gamma
}

expr16 = sympy.sympify(formulas[16], locals=names)
f16 = sympy.lambdify(['speed', 'angle'], expr16, ['numpy', 'scipy'])

sympy can also write the formula in C, JavaScript, LaTeX and other languages:

print(sympy.ccode(expr))
# 0.10197858992410698*pow(speed, 2)*sin(1.9991926274176621*angle)

That formula runs on a microcontroller, in a spreadsheet, or in a simulation, with no model file and no Python.

Why a formula and not gradient boosting?

On data like the training data, a good black-box model does fine too. Gradient boosting from scikit-learn, trained on the same rows, scores R² = 0.988 on the test set, against 0.9998 for the formula. The difference shows up outside the range of the training data. Our throws had speeds from 5 to 30 m/s; here is what the two models predict for faster ones (speed in m/s, angle in radians, distances in meters):

from sklearn.ensemble import GradientBoostingRegressor

boost = GradientBoostingRegressor(random_state=0).fit(train[['speed', 'angle']], train.distance)
far = pd.DataFrame({
    'speed': [40.0, 50.0, 60.0],
    'angle': [0.6, 0.8, 1.0]
})
print(f(far.speed, far.angle), boost.predict(far))

Tree-based models can't extrapolate: outside the training data they predict about the largest values they saw (the longest training throw was 87 m), so the boosting predictions stay flat. The formula found the mechanism, so it keeps working. That is the case for symbolic regression: when the data comes from a law, finding the law beats approximating it. When it doesn't, the formula will be an approximation too, and the test score tells you how good it is.

Useful options

start_process accepts all of TuringBot's search settings. The ones you will use most:

• allowed_functions: the building blocks, for example '+ * / sqrt exp'. Restricting them to what you expect gives simpler, more plausible formulas.
• custom_formula: search for a formula of a given shape. With 'distance = f(angle)*speed*speed', the search only has to find how the distance depends on the angle; on this data it finds 0.102*sin(2.0*angle)*speed*speed in a few seconds (sim.functions shows the whole formula, with f filled in).
• search_metric: the error to minimize, RMS error (4) by default. Mean relative error (1) suits targets that span orders of magnitude; there are also classification metrics.
• train_test_split: TuringBot holds out part of the data (for example, 80 for 80/20) and reports a test error for each formula, as a third value in each row of sim.functions, to spot overfitting while it searches (so loop over complexity, error, test_error, formula).
• integer_constants=True: only integer constants, for formulas like x**2 + 3*y.
• threads: CPU threads; all of them by default. The free version uses up to 2.

To watch a search as it improves, call refresh_functions() in a loop instead of waiting, and stop it with terminate_process() when the formulas stop getting better. The full list of options is in the documentation.

If you prefer the scikit-learn interface

The library also has a scikit-learn estimator, TuringBotRegressor, for pipelines and tools like cross_val_score. It runs the same search and evaluates the chosen formula in predict:

from turingbot import TuringBotRegressor

model = TuringBotRegressor(max_time=30)
model.fit(train[['speed', 'angle']], train.distance)
print(model.sympy())  # The formula it chose
print(model.equations_[['complexity', 'loss', 'equation']])  # All of them
print(model.predict(test[['speed', 'angle']])[:3])

By default it picks the formula with the largest relative drop in error per unit of complexity, among those within 1.5 times the lowest error; on this data that is the same complexity 13 formula. It takes most of the search options above as parameters; custom_formula and a few others are only in start_process.

Other Python libraries for symbolic regression

Two open-source libraries are widely used: PySR and gplearn. Both have scikit-learn estimators (PySRRegressor and SymbolicRegressor) and use genetic programming, while TuringBot uses simulated annealing. Our guide to symbolic regression explains the methods and compares the tools.

Frequently asked questions

Does it work with NumPy arrays?

Yes. Pass the array as input_file, with the target in the last column, and its column names with column_names=['speed', 'angle', 'distance']; without them, the formulas use default names such as col1 and col2.

Does it work in Jupyter notebooks?

Yes. The search runs in a separate process, so the notebook stays usable while it runs; call refresh_functions() in a cell whenever you want to see the latest formulas.

What are the free version's limits?

200 rows, 5 columns and 2 CPU threads. The example above fits within them. A license removes the limits.