Statistical inference on a very small datasets

Question

I have been working with machine learning for about a year now, but mostly with large datasets. However, I am currently working on a problem with a very small dataset. Here is my problem: I am creating a rocket fuel with 4 ingredients, x1, x2, x3, x4, and I want to maximize reaction strength, y. I have already mixed them in the arrangements below to get the corresponding values.

1. (0.9)x1 + (0.0)x2 + (0.1)x3 + (0)x4 = 16.5
2. (0.0)x1 + (0.9)x2 + (0.1)x3 + (0)x4 = 8.6
3. (.45)x1 + (.45)x2 + (0.0)x3 + (0.1)x4 = 12.6
4. (0.6)x1 + (0.3)x2 + (0.05)x3 + (.05)x4 = 18.9

5. (0.3)x1 + (0.9)x2 + (0.05)x3 + (.05)x4 = 9.8

   My next question is, how should I design my next few mixtures to maximize the reaction strength? Can you suggest any algorithms or statistical frameworks to get me started? Much appreciated.

Answer 1

There are two separate issues:

1. Sampling - Picking the optional ingredient level for next experiment to run. Given you have only have 4 explanatory variables, just plot them. Either all pairwise or a couple of 3d charts. With the outcome variable on the y or z axis. You'll then see the trend in the data. You can decide to get more data for interpolation (between the data points you already have) or extrapolation (data outside of the current range). There are frameworks, such as Bayesian Optimization, but that is too much work given the small dimensionality.

2. Inference - Predicting performance for new data. Given the data you have seen thus far (sample data), estimate parameters. In your example that would the estimating the contribution of each of the 4 ingredients, either individually or interaction. Those parameters could be scalar coefficients or distributions.

Answer 2

This is a perfect problem for active learning. Methods based on Bayesian Optimization are particularly powerful for optimizing black-box functions which are expensive to evaluate (i.e. running an experiment in the lab). There are a few BO packages out there which may be of interest, Martin Kraisser's blog has a nice overview.

I noticed that the features in your last experiment don't add up to 1 which I am assuming was a typo. For the demo I changed that entry to x2 = 0.6.

Here is a sample I threw together in python using GPyOpt, a Gaussian Process based package:

import numpy as np
import GPyOpt

x_init = np.array([[0.9,0.0,0.1,0.0],
                   [0.0,0.9,0.1,0.0],
                   [0.45,0.45,0.0,0.1],
                   [0.6,0.3,0.05,0.05],
                   [0.3,0.6,0.05,0.05]])

y_init = np.array([[16.5],[8.6],[12.6],[18.9],[9.8]])*(-1)

domain = [{'name': 'x1', 'type': 'continuous', 'domain': (0,1.0)},
        {'name': 'x2', 'type': 'continuous', 'domain': (0,1.0)},
        {'name': 'x3', 'type': 'continuous', 'domain': (0,1.0)},
        {'name': 'x4', 'type': 'continuous', 'domain': (0,1.0)}        
        ]

constraints = [
        {'name':'const_1', 'constraint': '(x[:,0] + x[:,1] + x[:,2] + x[:,3]) - 1 - 0.001'},
        {'name':'const_2', 'constraint': '1 - (x[:,0] + x[:,1] + x[:,2] + x[:,3]) - 0.001'}
        ]

bo_step = GPyOpt.methods.BayesianOptimization(
        f = None,
        domain = domain,
        constraints = constraints,
        X = x_init,
        Y = y_init
        )

x_next = bo_step.suggest_next_locations()

print(x_next)
print(np.sum(x_next))

Notes:

(1) GPyOpt only accepts constraints in a certain form, that's why there are 2 which constrain y on the interval [0.999,1.001].

(2) GPyOpt is set up to minimize functions. Therefore I multiplied your target values by -1.

This example suggests that your next experiment should be run at:

x1 = 0.69 x2 = 0.21 x3 = 0.08 x4 = 0.02

BO algorithms can be tuned to give different results based on your preferences for exploiting existing information versus exploring new areas of the space. I'm not sure what GPyOpts standard settings are so If you are interested it could be worth looking at the documentation.