Bayesian Machine Scientist¶

Example¶

In [ ]:

# Uncomment the following line when running on Google Colab
# !pip install autora

# Uncomment the following line when running on Google Colab # !pip install autora

Let's generate a simple data set with two features $x_1, x_2 \in [0, 1]$ and a target $y$. We will use the following generative model: $y = 2 x_1 - e^{(5 x_2)}$

In [7]:

import numpy as np

x_1 = np.linspace(0, 1, num=10)
x_2 = np.linspace(0, 1, num=10)
X = np.array(np.meshgrid(x_1, x_2)).T.reshape(-1,2)

y = 2 * X[:,0] + np.exp(5 * X[:,1])

import numpy as np x_1 = np.linspace(0, 1, num=10) x_2 = np.linspace(0, 1, num=10) X = np.array(np.meshgrid(x_1, x_2)).T.reshape(-1,2) y = 2 * X[:,0] + np.exp(5 * X[:,1])

Now let us choose a prior over the primitives. In this case, we will use priors determined by Guimerà et al (2020).

In [8]:

prior = "Guimera2020"

prior = "Guimera2020"

Set up the BMS Regressor¶

We will use the BMS Regressor to predict the outcomes. There are a number of parameters that determine how the architecture search is performed. The most important ones are listed below:

• epochs: The number of epochs to run BMS. This corresponds to the total number of equation mutations - one mcmc step for each parallel-tempered equation and one tree swap between a pair of parallel-tempered equations.
• prior_par: A dictionary of priors for each operation. The keys correspond to operations and the respective values correspond to prior probabilities of those operations. The model comes with a default.
• ts: A list of temperature values. The machine scientist creates an equation tree for each of these values. Higher temperature trees are harder to fit, and thus they help prevent overfitting of the model.

Let's use the same priors over primitives that we specified on the previous page as well as an illustrative set of temperatures to set up the BMS regressor with default parameters.

In [9]:

from autora.skl.bms import BMSRegressor

temperatures = [1.0] + [1.04**k for k in range(1, 20)]

primitives = {
    "Psychology": {
        "addition": 5.8,
        "subtraction": 4.3,
        "multiplication": 5.0,
        "division": 5.5,
    }
}

bms_estimator = BMSRegressor(
    epochs=1500,
    prior_par=primitives,
    ts=temperatures,
)

from autora.skl.bms import BMSRegressor temperatures = [1.0] + [1.04**k for k in range(1, 20)] primitives = { "Psychology": { "addition": 5.8, "subtraction": 4.3, "multiplication": 5.0, "division": 5.5, } } bms_estimator = BMSRegressor( epochs=1500, prior_par=primitives, ts=temperatures, )

Now we have everything to fit and verify the model.

In [10]:

bms_estimator.fit(X,y)
bms_estimator.predict(X)

bms_estimator.fit(X,y) bms_estimator.predict(X)

INFO:autora.skl.bms:BMS fitting started
  0%|          | 0/1500 [00:00<?, ?it/s]

---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
Cell In[10], line 1
----> 1 bms_estimator.fit(X,y)
      2 bms_estimator.predict(X)

File ~/Developer/autora/autora/skl/bms.py:133, in BMSRegressor.fit(self, X, y, num_param, root, custom_ops, seed)
    120     self.add_primitive(root)
    121 self.pms = Parallel(
    122     Ts=self.ts,
    123     variables=self.variables,
   (...)
    131     seed=seed,
    132 )
--> 133 self.model_, self.loss_, self.cache_ = utils.run(self.pms, self.epochs)
    134 self.models_ = list(self.pms.trees.values())
    136 _logger.info("BMS fitting finished")

File ~/Developer/autora/autora/theorist/bms/utils.py:35, in run(pms, num_steps, thinning)
     33 desc_len, model, model_len = [], pms.t1, np.inf
     34 for n in tqdm(range(num_steps)):
---> 35     pms.mcmc_step()
     36     pms.tree_swap()
     37     if num_steps % thinning == 0:  # sample less often if we thin more

File ~/Developer/autora/autora/theorist/bms/parallel.py:102, in Parallel.mcmc_step(self, verbose, p_rr, p_long)
     99     p_rr = 0.0
    100 for T, tree in list(self.trees.items()):
    101     # MCMC step
--> 102     tree.mcmc_step(verbose=verbose, p_rr=p_rr, p_long=p_long)
    103 self.t1 = self.trees["1.0"]

File ~/Developer/autora/autora/theorist/bms/mcmc.py:1160, in Tree.mcmc_step(self, verbose, p_rr, p_long)
   1157 else:
   1158     # Try to replace the root
   1159     newrr = choice(self.rr_space)
-> 1160     dE, dEB, dEP, par_valuesNew = self.dE_rr(rr=newrr, verbose=verbose)
   1161     if self.num_rr > 0 and -dEB / self.BT - dEP / self.PT > 0:
   1162         paccept = 1.0

File ~/Developer/autora/autora/theorist/bms/mcmc.py:1093, in Tree.dE_rr(self, rr, verbose)
   1090 self.par_values = old_par_values
   1092 # Prior: change due to the numbers of each operation
-> 1093 dEP += self.prior_par["Nopi_%s" % rr[0]]
   1094 try:
   1095     dEP += self.prior_par["Nopi2_%s" % rr[0]] * (
   1096         (self.nops[rr[0]] + 1) ** 2 - (self.nops[rr[0]]) ** 2
   1097     )

KeyError: 'Nopi_*'

Troubleshooting¶

We can troubleshoot the model by playing with a few parameters:

• Increasing the number of epochs. The original paper recommends 1500-3000 epochs for reliable fitting. The default is set to 1500.
• Using custom priors that are more relevant to the data. The default priors are over equations nonspecific to any particular scientific domain.
• Increasing the range of temperature values to escape local minima.
• Reducing the differences between parallel temperatures to escape local minima.