---
abstract: 'We introduce a method for non-uniform random number generation based on sampling a physical process in a controlled environment. We demonstrate one proof-of-concept implementation of the method that reduces the error of Monte Carlo integration of a univariate Gaussian by $1068 \times$ while doubling the speed of the Monte Carlo simulation. We show that the supply voltage and temperature of the physical process must be controlled to prevent the mean and standard deviation of the random number generator from drifting.'
author:
- 'James Timothy Meech, Phillip Stanley-Marbell,  [^1]'
bibliography:
- 'ms.bib'
title: Efficient Programmable Random Variate Generation Accelerator from Sensor Noise
---

Sensor, Noise, Bayesian, Inference, Non-uniform, Random.

Introduction {#section:introduction}
============

software-based methods of non-uniform random variate generation are slow and inefficient [@CCPUGPUFPGAMPPARNG][@NRETBDPC][@NURNGTPLA][@EHGRVAD]. We present a programmable system capable of generating Gaussian random variates by extracting the noise properties of a MEMS sensor. We demonstrate its principle and application. Sampling a random physical process with a given distribution provides a continuous random variable with a theoretically unlimited sample rate. In a hardware implementation, the sample rate of the analog to digital convertor limits random number generation rate.

Generating Uniform Random Variates Is Easy {#easy}
------------------------------------------

Uniform random numbers are generated for cryptography [@RRS]. Gaussian random variate generation is typically an order of magnitude slower and less efficient than uniform random variate generation [@CCPUGPUFPGAMPPARNG]. Table \[Table:gaussian\] compares several state-of-the-art methods of generating non-uniform random variates with a Gaussian distribution. We propose a method superior to all of the state-of-the-art methods in terms of sample rate and efficiency, consisting of a physical noise source and an analog to digital convertor.

  Device    Sample rate (Gsample/s)   Efficiency (Msample/J)  
  -------- ------------------------- ------------------------ --
  CPU                0.89                      3.17           
  GPU                12.90                    107.52          
  MPPA               0.86                     61.48           
  FPGA               12.10                    403.20          
  ADC                14.00                    645.10          

  : Comparison of state-of-the-art methods for 16-bit Gaussian random variate generation [@CCPUGPUFPGAMPPARNG]. The final row of the table represents seven dual-channel 1 Gsample/s ADS54J60s sampling each axis of five ADXL335s. Assuming that the samples move into the CPU memory in negligible time.[]{data-label="Table:gaussian"}

Generating Non-Uniform Random Variates Is Harder {#hard}
------------------------------------------------

The inversion method and accept reject method are used in software for generating samples from non-uniform random variates [@NURVG]. Let $U$ and $X$ be uniform and non-uniform random variates respectively and $F^{-1}$ the analytical closed-form solution for the inverse cumulative distribution function [@NURVG]. Algorithm \[Algorithm1\] shows the inversion method. The inversion method requires that the inverse cumulative distribution function has an analytical closed-form solution. The Gaussian distribution has no analytical closed-form solution for its inverse cumulative distribution function so cannot be used with the inversion method [@EHGRVAD]. The accept reject method must be used instead. Let $U$ and $X$ be independent and $g$ the density on $\mathbb{R}^d$ of $X$. Where $\mathbb{R}^d$ is $d$-dimensional Euclidean space. Algorithm \[Algorithm2\] shows the accept reject method. The accept reject method requires more math operations than the inversion method to transform a sample [@NURVG]. This causes it to take more clock cycles to compute and therefore, more time to transform each sample. When using the accept reject method, samples deviating from the desired distribution are rejected [@NURVG].

Generate uniform \[0,1\] random variate $U$ RETURN $X \leftarrow F^{-1}(U)$

initialization $T=1$ initialization $U=1$

Uses Of Non-Uniform Random Variates
-----------------------------------

Non-uniform random variate generators are fundamental to applications employing Monte Carlo Methods [@REGFPUED][@MCMIT][@PCRAPA], population balance modelling of the crystallisation process [@CQMMMCAESPBMDMCSDCCP], generating phylogenetic trees [@MLPRMA], drug discovery [@MCSPICACS], ray tracing [@SSCG], communication channel emulation [@EFPGAIGNGCCE], financial computing [@RAMCBFS], local area networks [@SLAN], modelling of manufacture systems [@SMMS] and measuring healthcare strategy effectiveness [@BBIGHSE].\
When conducting Bayesian inference in probabilistic machine learning, we must evaluate Bayes’ theorem to calculate the probability that a belief $B$ is true given new data $D$, we denote this as $P(B|D)$. To do this we need the probability that the belief is true regardless of our data $P(B)$, the probability that the data is true given the belief $P(D|B)$ and the probability that the data is true regardless of the belief $P(D)$. We will refer to $P(B)$ as the prior, $P(D|B)$ as the likelihood, $P(D)$ as the marginal likelihood and $P(B|D)$ as the posterior [@lambert2018student].

$$P(B|D) = \frac{P(D|B)P(B)}{P(D)}.$$

We calculate the marginal likelihood by integrating the joint density $P(B,D)$ [@lambert2018student]. In practice, the analytical calculation of the marginal likelihood is impossible for all but the simplest joint distributions [@lambert2018student]. We instead sample from the joint distribution $P(B,D)$ to obtain summary statistics that we can use to describe it [@lambert2018student]. These random samples from bespoke probability distributions must be produced using a non-uniform random variate generator.\

Contributions
-------------

1.  The idea that physical noise sources such as MEMS sensors can be used as non-uniform random variate generators (Section \[section:introduction\]).

2.  Estimation of performance increase and error reduction achieved by using such a generator for Monte Carlo integration (Section \[section:MotivatingExample\]).

3.  Investigation of the impact of temperature and supply voltage on the noise distribution obtained from a commercial MEMS sensor (Section \[section:ThermalChamberExperiment\]).

Motivating Example {#section:MotivatingExample}
==================

We performed Monte Carlo integration of a Gaussian with a mean of 0 and standard deviation of 1 using samples from a Gaussian with the same parameters. We ran the experiment with a Gaussian generated by the C++ random library and repeated it with the random number generation rate adjusted to match that of our proposed non-uniform random variate generator based on physical noise sources. We performed the same integration using samples from a uniform distribution with a minimum of -3 and a maximum of 3. Let $e$ be the error of the integration, $t$ be the time taken by the integration, $N$ be the number of random samples and $D$ be the distribution (either uniform or Gaussian). Let samps be the array of random numbers, $A$ be the area, $b$ be the rectangle base length, $h$ be the rectangle height and $f$ the probability density function of the Gaussian for integration. Algorithm \[Algorithm3\] shows the integration scheme that we used. We repeated each integration process 1000 times and calculated the average. Figure \[fig:MonteError\] shows that the proposed hardware random number generator outperforms the C++ uniform random number generator after 1000 samples and reduces the error by $1068 \times$ after 1 million samples. The proposed hardware random number generator always performs the task at least twice as fast as the C++ Gaussian random number generator.

Timer start Generate $N$ random samples from distribution $D$ Sort $N$ random samples e = abs(1-A) Timer stop t = stop - start RETURN e,t

![The error of Monte Carlo integration plateaus for uniform random numbers but decreases with a power relation for Gaussian random numbers.[]{data-label="fig:MonteError"}](Figures/FinalMonteError.pdf){width="45.00000%"}

Methodology {#section:methodology}
===========

A non-uniform random variate generator based on physical noise sources must have negligible drift of the mean and standard deviation over time. Drift would cause errors in calculations using the output of the random number generator. Any environmental parameter that causes non-negligible drift must be physically controlled, for example using a voltage regulator to control voltage and a Peltier device to control temperature. We investigated the dependence of 100,000 sample distributions on voltage and temperature. We sampled the z-axis of a MEMS accelerometer (we used the accelerometer in the Bosch BMX055) to obtain the distributions.

Temperature-Controlled Experiments {#section:ThermalChamberExperiment}
----------------------------------

Figure \[Figure:CircuitDiagram\] shows the experimental setup. We placed the microcontroller, accelerometer, tilt and rotate stage and vibration isolation platform inside a Binder MK56 thermal chamber. We connected a microcontroller to the sensor via I2C for a 1154 Hz sample rate. We used a Keithly 2450 source measurement unit to power the sensor and measure the current drawn. We set the chamber temperature to 25 $^\circ$C and allowed 30 minutes for the temperature of the sensor to equilibrate whilst constantly sampling z-axis acceleration values from it. We then sampled 100,000 values from the BMX055 sensor at a 3.6 V supply voltage. We repeated this for all the voltages in the range of 3.6 to 1.4 V with a 0.2 V decrement. We then repeated this process for the temperature from 25 down to -5 $^\circ$C with a decrement of 5 $^\circ$C.

![Experimental setup, powering the sensor and the I2C interface separately allowed more accurate measurement of the power consumed by the sensor.[]{data-label="Figure:CircuitDiagram"}](Figures/CircuitDiagram.pdf){width="45.00000%"}

Quantisation Investigation
--------------------------

We investigated the effect of quantisation on the Kullback–Leibler (KL) divergence between a discrete distribution and its ideal fitted curve. We used the MATLAB `normrnd` function to generate 100,000 values from a Gaussian distribution with the same mean and standard deviation as the BMX055 z-axis at 2.6 V and 10 $^\circ$C. We then discretized the values into various numbers of bins, fitted a Gaussian distribution to them and calculated the KL divergence between the fitted distribution and the actual distribution.

Results and Discussion
======================

We calculated the KL divergence between two discrete distributions using the following equation where $P$ and $Q$ are discrete probability distributions, $x$ is a given sample value and $\chi$ the sample space [@KLDIV]:

$$D_{KL}(P||Q) = - \sum_{x \in \chi} P(x) \log \bigg(\frac{Q(x)}{P(x)} \bigg).$$

Figure \[Figure:Gauss\] shows the BMX055 z-axis acceleration distribution at 2.6 V and 10 $^\circ$C. We found that the KL divergence between the distribution from the BMX055 z-axis and its fitted Gaussian (0.00263) was more than an order of magnitude smaller than the equivalent result for a MATLAB generated distribution (0.0392). We rounded the MATLAB generated floats for comparison with the BMX055 generated integers. The KL divergence between a MATLAB generated uniform distribution and its fitted Gaussian is 0.116 for reference. We averaged the KL divergence calculations over 100 distributions to account for the random variations in the measurement. The numbers generated by the sensor are closer to an ideal Gaussian distribution than those generated by MATLAB.

![Distribution of BMX055 accelerometer z-axis noise and closest fitted distribution. The KL divergence or difference between the data and the fitted distribution is 0.00263 demonstrating that the accelerometer noise closely matches a Gaussian distribution.[]{data-label="Figure:Gauss"}](Figures/Distribution.pdf){width="45.00000%"}

Figures \[Figure:trendmean\] and \[Figure:trendsd\] show how voltage and temperature effect the mean and standard deviation of the BMX055 z-axis acceleration measurement. Temperature has a greater effect on mean and standard deviation than voltage. Both voltage and temperature have a greater effect upon the standard deviation than the mean. A random number generator based on this phenomenon should control both the temperature and the supply voltage of the sensor.

![The mean of BMX055 z-axis distributions increases with increasing temperature and supply power. Normalised to the maximum observed mean.[]{data-label="Figure:trendmean"}](Figures/MeanPercent.pdf){width="45.00000%"}

![The standard deviation of BMX055 z-axis distributions decreases with increasing temperature and supply power. Normalised to the maximum observed standard deviation.[]{data-label="Figure:trendsd"}](Figures/SDPercent.pdf){width="45.00000%"}

Figure \[Figure:discret\] shows the effect of increasing the bin size on the divergence between a distribution of 100,000 values and its ideal fitted distribution. This shows that increased quantization increases the difference between a distribution and its fitted Gaussian.

![The effect of increasing discretisation on the KL divergence whilst keeping the total number of samples constant at 100,000. Increased quantisation increases KL divergence.[]{data-label="Figure:discret"}](Figures/KLDIV.pdf){width="45.00000%"}

Gaussian To Gaussian Transform
==============================

A univariate Gaussian can be transformed to any other univariate Gaussian with one multiplication and one addition. This is significantly less computation than the accept reject method which requires at least 10 operations per each accept reject test: an exponential, square, square root, subtraction, comparison and five divisions / multiplications. The accept reject method may furthermore need to repeat this set of operations numerous times for each random variate. Two uniform random numbers must be generated for one accept reject method sample. Rapid addition and multiplication can be achieved using fast adders and multipliers implemented on an FPGA, Figure \[Figure:GausstoGauss\] shows how this could be achieved. The CPU requests a distribution by specifying parameters to the transform circuitry which proceeds to transform the input Gaussian to fit the requested output distribution. The transform circuitry then stores the values in a small high-speed cache that the CPU can read from. The asynchronous transform circuitry on the FPGA will be free to execute transformations at a rate not bound by the global clock frequency. Offloading the transformation to dedicated hardware leaves the processor free to execute other instructions which will improve performance.

![Gaussian to Gaussian transform implementation. The CPU can request distributions with certain parameters and they are presented to it in a small fast cache.[]{data-label="Figure:GausstoGauss"}](Figures/GausstoGauss.pdf){width="45.00000%"}

Conclusion {#section:summary}
==========

Sensors are a feasible source of non-uniform random variates at a higher sample rate and with greater efficiency than all of the state-of-the-art methods used to generate them. The mean and standard deviation of the noise produced by the z-axis of a commercial accelerometer depend upon the temperature of the environment and the supply voltage. The parameters of the distribution drift with temperature and supply voltage so both must be kept constant to avoid this. Quantizing a Gaussian distribution increases the KL divergence between it and a fitted Gaussian. This type of non-uniform random number generator can reduce the error of Monte Carlo integration of a univariate Gaussian by a factor of $1068 \times$ whilst doubling the speed.

[^1]: James Timothy Meech and Phillip Stanley-Marbell are with the Department of Electrical Engineering, University of Cambridge, Cambridge, CB3 0FA UK e-mail: jtm45@cam.ac.uk.
