The application of adaptive measurement methods have a long tradition in psychophysics that allow to optimize the measurement process by placing stimuli such that the expected information gain is maximized. This page is about the method described in (Poppe, Benner, & Elze, 2012).

Examples

The source package contains several examples contained in the data/ subdirectory including the computation of binning posteriors, adaptive sampling simulations and experiments, and spike-train analysis. We discuss two of them in the following:

Python Interface

The software package contains several examples that can be tested directly after the software has been installed. The package provides the following commands:

  • bayesian-binning to compute posterior quantities which can be either saved to a configure file or plotted with matplotlib,
  • adaptive-sampling for adaptive sampling in real experiments or simulations.

Example 1 (bayesian-binning)

Assume that we have the results of an experiment with 16 stimuli and two possible responses, e.g. success and failure. To visualize the resulting posterior we save the experimental outcome in a file called data1.cfg:

[Counts]
counts:
  10 11 10 10 12 10 11 20 21 19 19 20 20 19 21 20
  90 89 90 90 88 90 89 80 79 81 81 80 80 81 79 80
parameters: data1-parameters.py
visualization: data1-visualization.py

The counts give the number of successes or failes for all 16 stimuli. The prior parameters are defined in the python script data1-parameters.py:

def parameters(K, L):
	alpha = generate_alpha(np.ones([K, L]))
	beta  = generate_beta([], L)
	gamma = generate_gamma(L)

	return alpha, beta, gamma

The functions generate_alpha, generate_beta, and generate_gamma are provided by the library and compute a default set of parameters. Here, L=16 is the number of stimuli and K=2 the number of responses, which the library determines from the counts matrix in the data1.cfg configuration file. To generate plots we specify the axis labels in data1-visualization.py:

def preplot(result, options):
	x = np.arange(0, len(result['moments'][0]), 1)
	return x, ``

def postplot(ax, p, result, options):
	[ax11, ax12, ax21, ax22] = ax 
	[ p11,  p12,  p21,  p22] = p
	if result['bprob'] and options['bprob']:
		ax12.set_ylabel(r'$P(\Rsh_i|E)$', font)
	ax11.set_ylabel(r'$P(s_i|E)$', font)
	ax11.set_xlabel(r'$x$',  font)
	ax21.set_ylabel(r'$P(M = m_B|E)$', font)
	ax21.set_xlabel(r'$m_B$',  font)

Now run

bayesian-binning -v -b -m data1.cfg

to plot the posterior including break-probabilities and posterior density. To save the plot to a file use

bayesian-binning -v -b -m --savefig=data1.png data1.cfg

The result should look as follows: [image1]: images/binning.png ![Bayesian binning marginal posterior plot.][image1]

The figure shows the marginal posterior plot. The upper plot shows the expectation (thick red line), standard deviation (thin black line), and marginal density (gray shadings). The green line shows the break probability. The lower plot shows the model posterior, i.e. the posterior probability of multibins with m bins

Example 2 (adaptive-sampling)

Assume we want to do an experiment with L=36 stimuli and K=2 possible responses, i.e. success and failure. We can simulate such an experiment with the adaptive-sampling with a given ground truth that determines the input-response relation. The choose the configuration

[Ground Truth]
gt: 0.983707 0.977895 0.970075 0.959606 0.945683 0.927331 0.903428 0.872769 0.834219 0.786948 0.730763 0.666468 0.59614 0.523201 0.452233 0.388603 0.338142 0.307012 0.301628 0.327954 0.38205 0.481978 0.593685 0.705735 0.901878 0.924494 0.944262 0.956063 0.975398 0.986731 0.993179 0.996743 0.998648 0.999621 1.00000
visualization: data2-visualization.py

which we save as data2.cfg. The gt option sets the ground truth, that is, for each of the 36 stimuli the numeric value gives the probability of response success. The data2-visualization.py is given by

def preplot(result, options):
	x = np.arange(1, 36, 1)
	return x, ""

def postplot(ax, p, result, options):
	import re
	[ax11, ax12, ax21, ax22, ax31, ax32] = ax
	[ p11,  p12,  p21,  p22,  p31,  p32] = p
	ax11.lines[2].set_linewidth(2.5)
	ax12.lines[0].set_linewidth(0.0)
	ax12.lines[0].set_linestyle('--')
	ax12.lines[0].set_marker('_')
	ax12.lines[0].set_markersize(12)
	ax12.lines[0].set_mew(2.5)
	ax11.collections[0].set_alpha(0.05)
	ax11.set_xlim(1,35)
	ax11.set_xlabel(r'$x$')
	ax11.set_ylabel(r'$p_{x,s}$')
	ax11.set_ylim(0,1)
	ax12.set_ylabel('Ground truth')
	ax21.set_xlim(1,35)
	ax21.set_ylim(0,30)
	ax21.set_ylabel('Counts')
	ax21.set_xlabel(r'$x$')
	ax21.legend([p22], ['$U(x)$'], loc=4)
	ax22.set_ylabel(r'$U(x)$')

To run the simulation we type

adaptive-sampling -n 200 -m --no-model-posterior data2.cfg

The result should look like:

Adaptive sampling posterior plot.

The figure shows the marginal posterior after 200 samples. The lower plot shows the number of samples for each stimulus.

GNU-R Interface

We provide two GNU-R packages, one for adaptive sampling and another one that provides a fast implementation of the proximal multibin summation (ProMBS) algorithm.

Example 1 (adaptive.sampling)

Assume that we have an experiment with L=6 stimuli and K=2 responses. A possible configuration could look like

require('adaptive.sampling')

L = 6 # number of stimuli
K = 2 # number of responses
counts_success <- c(2,3,2,4,7,7)
counts_failure <- c(8,7,7,6,3,2)
counts <- count.statistic(t(matrix(c(counts_success, counts_failure), L)))
alpha_success  <- c(1,1,1,1,1,1)
alpha_failure  <- c(1,1,1,1,1,1)
alpha  <- default.alpha(t(matrix(c(alpha_success, alpha_failure), L)))
beta   <- default.beta(L)
gamma  <- default.gamma(L)

where counts_success and counts_failure contain the counts of previous measurements. To obtain the sampling utility we call the library with

> sampling.utility(counts, alpha, beta, gamma)
$utility
[1] 0.02317324 0.01987610 0.02100055 0.02301429 0.02601794 0.02921347

$expectation
    [,1]      [,2]      [,3]      [,4]      [,5]      [,6]
[1,] 0.2721838 0.2960228 0.2971329 0.3972445 0.6612351 0.7005378
[2,] 0.7278162 0.7039772 0.7028671 0.6027555 0.3387649 0.2994622

attr(,"class")
[1] "sampling.utility"

The field $utility gives the expected utility for each of the 6 stimuli, where we select the maximum to determine the next stimulus that we present during the experiment. The $expectation field can be ignored. It gives the posterior probability for each stimulus-response pair given our current observations. To visualize the result we use

result <- binning.posterior(counts, alpha, beta, gamma)
plot(result)

Example 2 (multibin.summation)

The multibin.summation package provides fast implementations of the prombs and extended prombs algorithm. A call to the prombs algorithm might look like

require(multibin.summation)

g <-  c(1, 2, 3, 4, 5)
f <-  c(1, 2, 3, 4, 5,
	0, 1, 2, 3, 4,
	0, 0, 1, 2, 3,
	0, 0, 0, 1, 2,
	0, 0, 0, 0, 1)

dim(f) <- c(5,5)
f      <- t(f)

exp(prombs(log(g), log(f)))

which results in

[5 40 63 32 5]

Matlab Interface

The Matlab package provides the following commands:

  • binningPosterior: compute marginal posterior quantities
  • samplingUtility: compute expected utilties for the optimization of the experimental process
  • prombs: call the multibin summation algorithm
  • prombsExtended: call the extended multibin summation algorithm

Example 1 (sampling)

Assume a simple experiment with L stimuli and K possible responses. The following code shows a possible call of samplingUtiltiy:

bare_counts = [ 29 6 5 6 7 9 16 22 29 20 12 9 10 18 9 4 3 3 3 6 10 16 19 21 24 17 22 11 11 7 6 6 5 6 27;
                 0 0 0 0 0 0 1 2 6 8 5 6 6 11 20 8 6 7 7 13 14 13 12 8 6 4 2 0 1 0 0 0 0 0 0];

% K: responses
% L: stimuli
[K, L]  = size(bare_counts);

counts  = count_statistic(bare_counts);
alpha   = default_alpha(ones(K, L));
beta    = default_beta(L);
gamma   = default_gamma(L);

options = default_options();

result  = samplingUtility(counts, alpha, beta, gamma, options);

bare_counts contains the experimental outcomes of previous measurements. The expected utility for each stimulus is then given by result.utility. To select the next stimulus for our experiment we would simply compute the argmax of this vector. The binning posterior can be computed and visualized with

result = binningPosterior(counts, alpha, beta, gamma, options)
plot_binning(result)

References

  1. Poppe, S., Benner, P., & Elze, T. (2012). A predictive approach to nonparametric inference for adaptive sequential sampling of psychophysical experiments. J Math Psychol, 56(3), 179–195.