<jupyter_start><jupyter_text>Information maximiser Using neural networks, sufficient statistics can be obtained from data by maximising the Fisher information.The neural network takes some data ${\bf d}$ and maps it to a compressed summary $\mathscr{f}:{\bf d}\to{\bf x}$ where ${\bf x}$ can have the same size as the dimensionality of the parameter space, rather than the data space.To train the neural network a batch of simulations ${\bf d}_{\sf sim}^{\sf fid}$ created at a fiducial parameter value $\boldsymbol{\theta}^{\rm fid}$ are compressed by the neural network to obtain ${\bf x}_{\sf sim}^{\sf fid}$. From this we can calculate the covariance ${\bf C_\mathscr{f}}$ of the compressed summaries. We also use simulations ${\bf d}_{\sf sim}^{\sf fid+}$ created above the fiducial parameter value $\boldsymbol{\theta}^{\sf fid+}$ and simulations ${\bf d}_{\sf sim}^{\sf fid-}$ created below the fiducial parameter value $\boldsymbol{\theta}^{\sf fid-}$. These are compressed using the network and used to find mean of the summaries $\partial\boldsymbol{\mu}_\mathscr{f}/\partial\theta_\alpha\equiv\boldsymbol{\mu}_\mathscr{f},_\alpha$ where the numerical derivative is $({\bf x}_{\sf sim}^{\sf fid+}-{\bf x}_{\sf sim}^{\sf fid-})/(\boldsymbol{\theta}^{\sf fid+}-\boldsymbol{\theta}^{\sf fid-})$. We then use ${\bf C}_\mathscr{f}$ and $\boldsymbol{\mu}_\mathscr{f},_\alpha$ to calculate the Fisher information$${\bf F}_{\alpha\beta} = \boldsymbol{\mu}_\mathscr{f},^T_{\alpha}{\bf C}^{-1}_\mathscr{f}\boldsymbol{\mu}_\mathscr{f},_{\beta}.$$We want to maximise the Fisher information so to train the network we use the loss function$$\Lambda = -\frac{1}{2}|{\bf F}_{\alpha\beta}|.$$ When using this code please cite arXiv:1802.03537.The code in the paper can be downloaded as v1 or v1.1 of the code kept on zenodo:[![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.1175196.svg)](https://doi.org/10.5281/zenodo.1175196)The code presented below is version two (and is much more powerful) and is under constant development. This code is run using>`python-3.6.6`>`jupyter-4.4.0`>`tensorflow-1.10.1`>`numpy-1.14.5`>`tqdm==4.25.0`>`sys (native)`Although these precise versions may not be necessary, I have put them here to avoid possible conflicts.   Load modules<jupyter_code>%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import tensorflow as tf
import IMNN<jupyter_output>/Users/charnock/.pyenv/versions/3.6.6/lib/python3.6/importlib/_bootstrap.py:219: RuntimeWarning: compiletime version 3.5 of module 'tensorflow.python.framework.fast_tensor_util' does not match runtime version 3.6
  return f(*args, **kwds)<jupyter_text>Generate dataIn this example we are going to use $n_{\bf d}=10$ data points of a 1D field of Gaussian noise with zero mean and unknown variance to see if the network can learn to summarise this variance.We start by defining a function to generate the data with the correct shape. The shape must be```data_shape = None + input shape```<jupyter_code>input_shape = [10]<jupyter_output><empty_output><jupyter_text>It is useful to define the generating function so that it only takes in the value of the parameter as its input since the function can then be used for ABC later.The data needs to be generated at a fiducial parameter value and at perturbed values just below and above the fiducial parameter for the numerical derivative.<jupyter_code>θ_fid = 1.
Δθpm = 0.1<jupyter_output><empty_output><jupyter_text>The data at the perturbed values should have the shape```perturbed_data_shape = None + number of parameters + input shape```The generating function is defined so that the fiducial parameter is passed as a list so that many simulations can be made at once. This is very useful for the ABC function later.<jupyter_code>def generate_data(θ, train = False):
    if train:
        return np.moveaxis(np.random.normal(0., np.sqrt(θ), [1] + input_shape + [len(θ)]), -1, 0)
    else:
        return np.moveaxis(np.random.normal(0., np.sqrt(θ), input_shape + [len(θ)]), -1, 0)<jupyter_output><empty_output><jupyter_text>Training dataEnough data needs to be made to approximate the covariance matrix of the output summaries. The number of simulations needed to approximate the covariance is `n_s`. If the data is particularly large then it might not be possible to pass all the data into active memory at once and so several the simulations can be split into batches.For example if we wanted to make 2000 simulations, but estimate the covariance using 1000 simulations at a timewe would set```pythonn_s = 1000n_train = 2```We're going to use 1000 simulations to approximate the covariance and use only 1 combination<jupyter_code>n_s = 1000
n_train = 1<jupyter_output><empty_output><jupyter_text>The training data can now be made<jupyter_code>t = generate_data(θ = [θ_fid for i in range(n_train * n_s)], train = False)<jupyter_output><empty_output><jupyter_text>By suppressing the sample variance between the simulations created at the lower and upper parameter values for the numerical derivative, far fewer simulations are needed. We choose to use 5% of the total number of simulations.<jupyter_code>derivative_fraction = 0.05
n_p = int(n_s * derivative_fraction)<jupyter_output><empty_output><jupyter_text>The sample variance can be supressed by choosing the same initial seed when creating the upper and lower simulations.<jupyter_code>seed = np.random.randint(1e6)
np.random.seed(seed)
t_m = generate_data(θ = [θ_fid - Δθpm for i in range(n_train * n_p)], train = True)
np.random.seed(seed)
t_p = generate_data(θ = [θ_fid + Δθpm for i in range(n_train * n_p)], train = True)
np.random.seed()<jupyter_output><empty_output><jupyter_text>We also need to get the denominator of the derivative which is given by the difference between the perturbed parameter values    $$\frac{\partial}{\partial\theta} = \frac{1}{2\Delta\theta_{\pm}}.$$This needs to be done for every parameter and kept in a numpy array of shape `[number of parameters]`.<jupyter_code>derivative_denominator = 1. / (2. * Δθpm)
der_den = np.array([derivative_denominator])<jupyter_output><empty_output><jupyter_text>The fiducial simulations and simulations for the derivative must be collected in a dictionary to be stored on the GPU or passed to the training function.<jupyter_code>data = {"x_central": t, "x_m": t_m, "x_p":t_p}<jupyter_output><empty_output><jupyter_text>Test dataWe should also make some test data, but here we will use only one combination `n_train = 1`. This needs adding to the dictionary<jupyter_code>tt = generate_data([θ_fid for i in range(n_s)])
seed = np.random.randint(1e6)
np.random.seed(seed)
tt_m = generate_data([θ_fid - Δθpm for i in range(n_p)], train = True)
np.random.seed(seed)
tt_p = generate_data([θ_fid + Δθpm for i in range(n_p)], train = True)
np.random.seed()
data["x_central_test"] = tt
data["x_m_test"] = tt_m
data["x_p_test"] = tt_p<jupyter_output><empty_output><jupyter_text>Data visualisationWe can plot the data to see what it looks like.<jupyter_code>fig, ax = plt.subplots(1, 1, figsize = (10, 6))
ax.plot(data['x_central'][np.random.randint(n_train * n_s)], label = "training data")
ax.plot(data['x_central_test'][np.random.randint(n_s)], label = "test data")
ax.legend(frameon = False)
ax.set_xlim([0, 9])
ax.set_xticks([])
ax.set_ylabel("Data amplitude");<jupyter_output><empty_output><jupyter_text>It is also very useful to plot the upper and lower derivatives to check that the sample variance is actually supressed since the network learns extremely slowly if this isn't done properly.<jupyter_code>fig, ax = plt.subplots(2, 1, figsize = (10, 10))
plt.subplots_adjust(hspace = 0)
training_index = np.random.randint(n_train * n_p)
ax[0].plot(data['x_m'][training_index, 0], label = "lower training data", color = 'C0')
ax[0].plot(data['x_p'][training_index, 0], label = "upper training data", color = 'C0', linestyle = 'dashed')
test_index = np.random.randint(n_p)
ax[0].plot(data['x_m_test'][test_index, 0], label = "lower test data", color = 'C1')
ax[0].plot(data['x_p_test'][test_index, 0], label = "upper test data", color = 'C1')
ax[0].legend(frameon = False)
ax[0].set_xlim([0, 9])
ax[0].set_xticks([])
ax[0].set_ylabel("Data amplitude")
ax[1].axhline(xmin = 0., xmax = 1., y = 0., linestyle = 'dashed', color = 'black')
ax[1].plot(data['x_p'][training_index, 0] - data['x_m'][training_index, 0], color = 'C0')
ax[1].plot(data['x_p_test'][test_index, 0] - data['x_m_test'][test_index, 0], color = 'C1')
ax[1].set_xlim([0, 9])
ax[1].set_xticks([])
ax[1].set_ylabel("Difference between derivative data amplitudes");<jupyter_output><empty_output><jupyter_text>Initiliase the neural network Define network parametersThe network works with a base set of parameters which are> `'verbose'` - `bool` - whether to print out diagnostics> `'number of simulations'` - `int` - the number of simulations to use in any one combination> `'differentiation fraction'` - `float` - a fraction of the simulations to use for the numerical derivative> `'fiducial θ'` - `array` - fiducial parameters in an array> `'derivative denominator'` - `array` - denominator of the numerical derivative for each parameter> `'number of summaries'` - `int` - number of summaries the network makes from the data> `'input shape'` - `int` or `list` - the number of inputs or the shape of the input if image-like input> `'preload data'` - `dict` or `None` - the training (and test) data to preload as a TensorFlow constant in a dictionary, no preloading is done if `None`> `'calculate MLE'` - `bool` - whether to calculate the maximum likelihood estimate> `'prebuild'` - `bool` - whether to get the network to build a network or to provided your own> `'save file'` - `string` - a file name to save the graph (not saved if wrong type or not given)```pythonparameters = {    'verbose': True,    'number of simulations': n_s,    'differentiation fraction': derivative_fraction,    'fiducial θ': np.array([θ_fid]),    'derivative denominator': der_den,    'number of summaries': 1,    'input shape': input_shape,    'preload data': data,    'calculate MLE': True,    'prebuild': False,}```The module can also build simple convolutional or dense networks (or a mixture of the two), which can be trigger by setting `'prebuild': True`. Several parameters are required to allow the network to build the network. These are> `'wv'` - `float` - the variance with which to initialise the weights. If this is 0 or less, the network will determine the weight variance according to He initialisation> `'bb'` - `float` - the constant value with which to initialise the biases> `'activation'` - `TensorFlow function` - a native tensorflow activation function> `'α'` - `float` or `int` - an additional parameter, if needed, for the tensorflow activation function> `'hidden layers'` - `list` - the architecture of the network. each element of the list is a hidden layer. A dense layer can be made using an integer where thet value indicates the number of neurons. A convolutional layer can be built by using a list where the first element is an integer where the number describes the number of filters, the second element is a list of the kernel size in the x and y directions, the third elemnet is a list of the strides in the x and y directions and the final element is string of 'SAME' or 'VALID' which describes the padding prescription.Here is an example of the IMNN which uses 1000 simulations per combination and 50 upper and 50 lower simulations per derivative for a model with one parameter where we require one summary which are preloaded as a TensorFlow constant where we want to have access to the precomputed maximum likelihood estimate. The module will build the network which takes in an input array of shape `[10]` and allows the network to decide the weight initialisation, initialises the biases at `bb = 0.1` and uses `tf.nn.leaky_relu` activation with a negative gradient parameter of `α = 0.01`. The network architecture has two fully connected hidden layers with 128 neurons in each layer. The graph is saved into a file in the `data` folder called `model.meta`.<jupyter_code>parameters = {
    'verbose': True,
    'number of simulations': n_s,
    'fiducial θ': np.array([θ_fid]),
    'derivative denominator': der_den,
    'differentiation fraction': derivative_fraction,
    'number of summaries': 1,
    'calculate MLE': True,
    'prebuild': True,
    'input shape': input_shape,
    'preload data': data,
    'save file': "data/model",
    'wv': 0.,
    'bb': 0.1,
    'activation': tf.nn.leaky_relu,
    'α': 0.01,
    'hidden layers': [128, 128]
}<jupyter_output><empty_output><jupyter_text>Self-defined networkA self defined network can be used instead of letting the module build the network for you. This function needs to take in two input tensors, the first is the shape of the input with `None` in the first axis and the second tensor is a tensorflow float (which will be the dropout). Since the weights need to be shared between several corresponding networks each set of trainable variables must be defined in its own scope. An example of the above network defined outside of the module is```pythondef network(input_tensor, dropout):    with tf.variable_scope('layer_1'):        weights = tf.get_variable("weights", [10, 128], initializer = tf.random_normal_initializer(0., 1.))        biases = tf.get_variable("biases", [128], initializer = tf.constant_initializer(0.1))    x = tf.matmul(input_tensor, weights)    x = tf.add(x, biases)    x = tf.nn.leaky_relu(x, 0.01)    x = tf.nn.dropout(x, dropout)    with tf.variable_scope('layer_2'):        weights = tf.get_variable("weights", [128, 128], initializer = tf.random_normal_initializer(0., 1.))        biases = tf.get_variable("biases", [128], initializer = tf.constant_initializer(0.1))    x = tf.matmul(x, weights)    x = tf.add(x, biases)    x = tf.nn.leaky_relu(x, 0.01)    x = tf.nn.dropout(x, dropout)    x = tf.reshape(x, (-1, 300))    with tf.variable_scope('layer_3'):        weights = tf.get_variable("weights", [128, 1], initializer = tf.random_normal_initializer(0., np.sqrt(2. / 300)))        biases = tf.get_variable("biases", [1], initializer = tf.constant_initializer(0.1))    x = tf.matmul(x, weights)    x = tf.add(x, biases)    x = tf.nn.leaky_relu(x, 0.01)    x = tf.nn.dropout(x, dropout)    return x```  Initialise the network<jupyter_code>n = IMNN.IMNN(parameters = parameters)<jupyter_output>saving model as data/model.meta
network architecture is [[10], 128, 128, 1].<jupyter_text>Build the networkTo build the network a learning rate, η, must be defined.<jupyter_code>η = 1e-3<jupyter_output><empty_output><jupyter_text>The `setup(η)` function initialises the input tensors, builds the network and defines the optimisation scheme. If a self-defined network function (`network(a, b)`) has been constructed this can be passed to the setup function```pythonn.setup(η = η, network = network)```<jupyter_code>n.setup(η = η)<jupyter_output>Tensor("x:0", shape=(?, 10), dtype=float32)
Tensor("IMNN/layer_1/dense_1/mul:0", shape=(?, 128), dtype=float32)
Tensor("IMNN/layer_2/dense_2/mul:0", shape=(?, 128), dtype=float32)
Tensor("IMNN/layer_3/LeakyRelu:0", shape=(?, 1), dtype=float32)
Tensor("output:0", shape=(?, 1), dtype=float32)
Tensor("GatherNd:0", shape=(1000, 10), dtype=float32)
Tensor("IMNN_1/layer_1/dense_1/mul:0", shape=(1000, 128), dtype=float32)
Tensor("IMNN_1/layer_2/dense_2/mul:0", shape=(1000, 128), dtype=float32)
Tensor("IMNN_1/layer_3/LeakyRelu:0", shape=(1000, 1), dtype=float32)
Tensor("Reshape:0", shape=(50, 10), dtype=float32)
Tensor("IMNN_1/layer_1_1/dense_1/mul:0", shape=(50, 128), dtype=float32)
Tensor("IMNN_1/layer_2_1/dense_2/mul:0", shape=(50, 128), dtype=float32)
Tensor("IMNN_1/layer_3_1/LeakyRelu:0", shape=(50, 1), dtype=float32)
Tensor("Reshape_1:0", shape=(50, 10), dtype=float32)
Tensor("IMNN_1/layer_1_2/dense_1/mul:0", shape=(50, 128), dtype=float32)
Tensor("IMNN_1/layer_2_2/dense_2/mul:0", shape=[...]<jupyter_text>Changing minimisation schemeBy default the optimation scheme is```pythonn.backpropagate = tf.train.GradientDescentOptimizer(η).minimize(n.Λ)```To use any other training scheme, such as the `Adam` optimiser, it is sufficient to run```pythonn.backpropagate = tf.train.AdamOptimizer(η, β1, β2, ε).minimize(n.Λ)``` after `setup(η)` to override the default minimisation routine. Note that testing with Adam optimiser has found it to be incredibly unstable. If you want to continue to use the default minimisation routine but want to change the learning rate without reinitialising you can run ```pythonn.training_scheme(η = new_η)```  Train the networkWith the data we can now train the network. The function simply takes the number of epochs, `num_epochs`, the fraction of neurons kept when using dropout `keep_rate`, and the denominator for the derivative calculated earlier, `der_den`.<jupyter_code>num_epochs = 500
keep_rate = 0.8<jupyter_output><empty_output><jupyter_text>If the data has not been preloaded as a TensorFlow constant then it can be passed to the train function```pythontrain_F, test_F = n.train(num_epochs = num_epochs, n_train = n_train, keep_rate = keep_rate, data = data)```We can run<jupyter_code>n.train(num_epochs = num_epochs, n_train = n_train, keep_rate = keep_rate, data = data, history = True)<jupyter_output>100%|██████████| 500/500 [00:13<00:00, 36.29it/s, detF=5.34, detF_test=4.2]<jupyter_text>The train function automatically collects a dictionary of history elements when `history = True`. When `history = False` the dictionary only contains> `'F'` - the Fisher information at each epoch> `'det(F)'` - the determinant of the Fisher informationWhen `history = True` then the dictionary also contains> `'Λ'` - the loss function of the training data> `'μ'` - the mean of the fiducial simulations> `'C'` - the covariance of the fiducial simulations> `'det(C)'` - the determinant of the fiducial simulations> `'dμdθ'` - the mean of the numerical derivative of the simulationsTest version of each of these quantities is also calculated if there is provided test data> `'test F'` - the Fisher information of the test data> `'det(test F)'` - the determinant of the Fisher information from the test data> `'Λ'` - the loss function of the test data> `'μ'` - the mean of the fiducial test simulations> `'C'` - the covariance of the fiducial test simulations> `'det(C)'` - the determinant of the fiducial test simulations> `'dμdθ'` - the mean of the numerical derivative of the test simulations<jupyter_code>fig, ax = plt.subplots(5, 1, sharex = True, figsize = (8, 14))
plt.subplots_adjust(hspace = 0)
end = len(n.history["det(F)"])
epochs = np.arange(end)
a, = ax[0].plot(epochs, n.history["det(F)"], label = 'Training data')
b, = ax[0].plot(epochs, n.history["det(test F)"], label = 'Test data')
ax[0].legend(frameon = False)
ax[0].set_ylabel(r'$|{\bf F}_{\alpha\beta}|$')
ax[1].plot(epochs, n.history["Λ"])
ax[1].plot(epochs, n.history["test Λ"])
ax[1].set_xlabel('Number of epochs')
ax[1].set_ylabel(r'$\Lambda$')
ax[1].set_xlim([0, len(epochs)]);
ax[2].plot(epochs, n.history["det(C)"])
ax[2].plot(epochs, n.history["det(test C)"])
ax[2].set_xlabel('Number of epochs')
ax[2].set_ylabel(r'$|{\bf C}|$')
ax[2].set_xlim([0, len(epochs)]);
ax[3].plot(epochs, np.array(n.history["dμdθ"]).reshape((np.prod(np.array(n.history["dμdθ"]).shape))))
ax[3].plot(epochs, np.array(n.history["test dμdθ"]).reshape((np.prod(np.array(n.history["test dμdθ"]).shape))))
ax[3].set_ylabel(r'$\partial\mu/\partial\theta$')
ax[3].set_xlabel('Number of epochs')
ax[3].set_xlim([0, len(epochs)])
ax[4].plot(epochs, np.array(n.history["μ"]).reshape((np.prod(np.array(n.history["μ"]).shape))))
ax[4].plot(epochs, np.array(n.history["test μ"]).reshape((np.prod(np.array(n.history["test μ"]).shape))))
ax[4].set_ylabel('μ')
ax[4].set_xlabel('Number of epochs')
ax[4].set_xlim([0, len(epochs)])
print()<jupyter_output><jupyter_text>We can see that the test loss deviates from the training loss. This is to be expected because there are will be a lot of correlation within a small training set which isn't in the test set. As long as the test loss doesn't start increasing then it is likely that the network is still working, with the maximum Fisher available from the network is the value obtained from the test set.  Resetting the networkIf you need to reset the weights and biases for any reason then you can call```pythonn.reinitialise_session()```  Saving the networkIf you don't initialise the network with a save name you can save the network as a `TensorFlow` `meta` graph. For example saving the model in the directory `/.data` called `saved_model.meta` can be done using the function```pythonn.save_network(file_name = "data/saved_model", first_time = True)```If `save file` is passed with a correct file name when initialising the module then the initialised network will be saved by```pythonn.begin_session()```and then saved at the end of training.  Loading the networkYou can load the network from a `TensorFlow` `meta` graph (from `/.data/saved_model.meta`) using the same parameter dictionary as used when originally training the network and then running```pythonn = IMNN.IMNN(parameters = parameters)n.restore_network()```Training can be continued after restoring the model - although the Adam optimiser might need to reacquaint itself.  Approximate Bayesian computationWe can now do ABC (or PMC-ABC) with our calculated summary. First we generate some simulated real data:<jupyter_code>real_data = generate_data(θ = [1.], train = False)<jupyter_output><empty_output><jupyter_text>We can plot this real data to see what it looks like.<jupyter_code>fig, ax = plt.subplots(1, 1, figsize = (10, 6))
ax.plot(real_data[0], label = "real data")
ax.legend(frameon = False)
ax.set_xlim([0, 9])
ax.set_xticks([])
ax.set_ylabel("Data amplitude");<jupyter_output><empty_output><jupyter_text>ABCWe now perform ABC by drawing 100000 random samples from the prior. We define the upper and lower bounds of a uniform prior to be 0 and 10. Only a uniform prior is implemented at the moment. From the samples we create simulations at each parameter value and feed each simulation through the network to get summaries. The summaries are compared to the summary of the real data to find the distances which can be used to accept or reject points.Because the simulations are created within the ABC function then the generation function must be passed. This is why the generator should be of the form defined above, which takes only a list of parameter values and returns a simulation at each parameter.If the data is not preloaded as a TensorFlow constant then the data can be passed to the function as```pythonθ, summary, s, ρ, F = n.ABC(real_data = real_data, prior = [0, 10], draws = 100000, generate_simulation = generate_data, at_once = True, data = data)```Here we can use<jupyter_code>θ, summary, s, ρ, F = n.ABC(real_data = real_data, prior = [0, 10], draws = 100000, generate_simulation = generate_data, at_once = True, data = data)<jupyter_output><empty_output><jupyter_text>If the simulations are going to be too large to make all at once the `at_once` option can be set to false which will create one simulation at a time.```pythonθ, summary, s, ρ, F = n.ABC(real_data = real_data, der_den = der_den, prior = [0, 10], draws = 100000, generate_simulation = generate_data, at_once = False)```  Accept or rejectIn ABC draws are accepted if the distance between the simulation summary and the simulation of the real data are "close", i.e. smaller than some ϵ value, which is chosen somewhat arbitrarily.<jupyter_code>ϵ = 1
accept_indices = np.argwhere(ρ < ϵ)[:, 0]
reject_indices = np.argwhere(ρ >= ϵ)[:, 0]<jupyter_output><empty_output><jupyter_text>Plot samplesWe can plot the output samples and the histogram of the accepted samples, which should peak around `θ = 1` (where we generated the real data). The monotonic function of all the output samples shows that the network has learned how to summarise the data.<jupyter_code>fig, ax = plt.subplots(2, 1, sharex = True, figsize = (10, 10))
plt.subplots_adjust(hspace = 0)
ax[0].scatter(θ[accept_indices] , s[accept_indices, 0], s = 1)
ax[0].scatter(θ[reject_indices], s[reject_indices, 0], s = 1, alpha = 0.1)
ax[0].plot([0, 10], [summary[0], summary[0]], color = 'black', linestyle = 'dashed')
ax[0].set_ylabel('Network output', labelpad = 0)
ax[0].set_xlim([0, 10])
ax[1].hist(θ[accept_indices], np.linspace(0, 10, 100), histtype = u'step', density = True, linewidth = 1.5, color = '#9467bd');
ax[1].set_xlabel('$\\theta$')
ax[1].set_ylabel('$\\mathcal{P}(\\theta|{\\bf d})$')
ax[1].set_yticks([]);<jupyter_output><empty_output><jupyter_text>There can be a lot of $\theta$ draws which are unconstrained by the network because no similar structures were seen in the data which is indicative of using too small of a small training set.  PMC-ABCPopulation Monte Carlo ABC is a way of reducing the number of draws by first sampling from a prior, accepting the closest 75% of the samples and weighting all the rest of the samples to create a new proposal distribution. The furthest 25% of the original samples are redrawn from the new proposal distribution. The furthest 25% of the simulation summaries are continually rejected and the proposal distribution updated until the number of draws needed accept all the 25% of the samples is much greater than this number of samples. This ratio is called the criterion. The inputs work in a very similar way to the `ABC` function above. If we want 1000 samples from the approximate distribution at the end of the PMC we need to set `num_keep = 1000`. The initial random draw (as in ABC above) initialises with `num_draws`, the larger this is the better proposal distribution will be on the first iteration.If the data is not preloaded as a TensorFlow constant then the data can be passed to the function as```pythonθ_, summary_, ρ_, s_, W, total_draws, F = n.PMC(real_data = real_data, prior = [0, 10], num_draws = 1000, num_keep = 1000, generate_simulation = generate_data, criterion = 0.1, data = data, at_once = True, samples = None)```Here we can use<jupyter_code>θ_, summary_, ρ_, s_, W, total_draws, F = n.PMC(real_data = real_data, prior = [0, 10], num_draws = 1000, num_keep = 1000, generate_simulation = generate_data, criterion = 0.1, at_once = True, samples = None, data = data)<jupyter_output>iteration = 14, current criterion = 0.09885330170027679, total draws = 36132, ϵ = 2.1938307881355286.<jupyter_text>If we want the PMC to continue for longer we can provide the output of PMC as an input as```pythonθ_, summary_, ρ_, s_, W, total_draws, F = n.PMC(real_data = real_data, der_den = der_den, prior = [0, 10], num_draws = 1000, num_keep = 1000, generate_simulation = generate_data, criterion = 0.01, data = data, at_once = True, samples = [θ_, summary_, ρ_, s_, W, total_draws, F])```Finally we can plot the accepted samples and plot their histogram.<jupyter_code>fig, ax = plt.subplots(2, 1, sharex = True, figsize = (10, 10))
plt.subplots_adjust(hspace = 0)
ax[0].scatter(θ_ , s_, s = 1)
ax[0].plot([0, 10], [summary[0], summary[0]], color = 'black', linestyle = 'dashed')
ax[0].set_ylabel('Network output', labelpad = 0)
ax[0].set_xlim([0, 10])
ax[0].set_ylim([np.min(s_), np.max(s_)])
ax[1].hist(θ_, np.linspace(0, 10, 100), histtype = u'step', density = True, linewidth = 1.5, color = '#9467bd');
ax[1].set_xlabel('θ')
ax[1].set_ylabel('$\\mathcal{P}(\\theta|{\\bf d})$')
ax[1].set_yticks([]);<jupyter_output><empty_output><jupyter_text>Maximum likelihood estimateWe can also calculate the first-order Gaussian approximation of the posterior on the parameter and find a maximum likelihood estimate.If the data is not preloaded as a TensorFlow constant then it can be passed using```pythonasymptotic_likelihood = n.asymptotic_likelihood(real_data = real_data, prior = np.linspace(0, 10, 1000).reshape((1, 1, 1000)), data = data)MLE = n.θ_MLE(real_data = real_data, data = data)```Here we will use<jupyter_code>asymptotic_likelihood = n.asymptotic_likelihood(real_data = real_data, prior = np.linspace(0, 10, 1000).reshape((1, 1, 1000)), data = data)
MLE = n.θ_MLE(real_data = real_data, data = data)
fig, ax = plt.subplots(1, 1, figsize = (10, 6))
ax.plot(np.linspace(0, 10, 1000), asymptotic_likelihood[0, 0], linewidth = 1.5)
ax.axvline(x = MLE[0, 0], ymin = 0., ymax = 1., linestyle = 'dashed', color = 'black')
ax.set_xlabel("θ")
ax.set_xlim([0, 10])
ax.set_ylabel('$\\mathcal{P}(\\theta|{\\bf d})$')
ax.set_yticks([]);<jupyter_output><empty_output><jupyter_text>Analytic posterior calculationWe know what the analytic posterior is for this model$$\mathcal{P}(\boldsymbol{\theta}|{\bf d}) = \frac{\displaystyle{\rm exp}\left[-\frac{1}{2\boldsymbol{\theta}}\sum_{i = 1}^{n_{\bf d}}d_i\right]}{(2\pi\boldsymbol{\theta})^{n_{\bf d}/2}}.$$We can there for plot this as a comparison.<jupyter_code>θ_grid = np.linspace(0.001, 10, 1000)
analytic_posterior = np.exp(-0.5 * np.sum(real_data**2.) / θ_grid) / np.sqrt(2. * np.pi * θ_grid)**10.
analytic_posterior = analytic_posterior / np.sum(analytic_posterior * (θ_grid[1] - θ_grid[0]))
fig, ax = plt.subplots(1, 1, figsize = (10, 6))
ax.plot(θ_grid, analytic_posterior, linewidth = 1.5, color = 'C1', label = "Analytic posterior")
ax.hist(θ_, np.linspace(0, 10, 100), histtype = u'step', density = True, linewidth = 1.5, color = '#9467bd', label = "PMC posterior");
ax.hist(θ[accept_indices], np.linspace(0, 10, 100), histtype = u'step', density = True, linewidth = 1.5, color = 'C2', label = "ABC posterior")
ax.plot(np.linspace(0, 10, 1000), asymptotic_likelihood[0, 0], color = 'C0', linewidth = 1.5, label = "Asymptotic Gaussian likelihood")
ax.axvline(x = MLE[0, 0], ymin = 0., ymax = 1., linestyle = 'dashed', color = 'black', label = "Maximum likelihood estimate")
ax.legend(frameon = False)
ax.set_xlim([0, 10])
ax.set_xlabel('θ')
ax.set_ylabel('$\\mathcal{P}(\\theta|{\\bf d})$')
ax.set_yticks([]);<jupyter_output><empty_output>