Random Normal Generation with Arduino Rev 1.0

normal_dist_gen1.0.c

/**
 * Abstract:
 * Using the ATMEGA328/Arduino in this post, we will use a naive method to generate normally distributed random numbers. 
 * Our first step was to implement density and cumulative normal distribution functions. After that, in order to obtain 
 * an approximate inverse cumulative density function, a lookup table is used and z-values are obtained using linear interpolation. 
 * Using the Shapiro-Wilki test, we validated the normality of a generated number list.
 */

#include "support.h"
/**
 * 
 */
float get_linspacing(float min, float max, int length_out){
  return fabs(max - min)/(length_out - 1);
}

/**
 * 
 */
float* linspace(float min, float max, int length_out){
  float* seq = (float*) malloc(length_out * sizeof(float));
  int counter = 0;
  float h = fabs(max - min)/(length_out - 1);
  for(float i = min; i <= max + h; i = i + h){
    //printf("%d, %g\n", counter, i); 
    seq[counter] = i;
    counter++;
  }
  return seq;
}



/**
 * 
 */
typedef struct linfind_t {
  float a0; float a1;
  int i; int j;
  int count;
} linfind_t;

/**
 * 
 */
linfind_t linfind_approx(float* seq, float linspacing, int seqsize, float findval){ 
  linfind_t results;
  results.count = 0;
  for(int i = 0; i < seqsize; i++){
    if(fabs(seq[i] - findval) <= 1e-9){
      //printf("Value %g in index %d closer to %g\n", seq[i], i, findval);
      results.count = 1;
      results.a0 =  0;
      results.i  = -1;                
      //
      results.a1 = seq[i];
      results.j  = i;               
      break;
    } 
    //
    if(fabs(seq[i] - findval) <= linspacing ){
      //printf("(%g, %g, %d)\n", findval, seq[i], i);
      if(results.count == 0){       
        results.a0 =  seq[i];
        results.i  = i;
      }
      //
      if(results.count >= 1){               
        results.a1 =  seq[i];
        results.j  = i;
      }
      results.count++;
    } 
  }
  //
  return results;
}


/**
 * 
 */
float integrate(float (*f)(float, float, float), 
                float a, float b, int n_steps, 
                float mu, float sigma){
  float h = fabs((b - a)/n_steps);
  //
  float fi = 0.0;
  fi = fi + (*f)(a, mu, sigma) + (*f)(b, mu, sigma);
  for (float i = a + h; i < b; i = i + h){
    fi = fi + 2*((*f)(i, mu, sigma));
  }
  //
  return (h/2.0)*fi;
}

/**
 * 
 */
float dnorm(float x, float mu, float sigma){
  return 1.0/(sigma*sqrt(2*PI))*exp(-0.5*pow((x-mu)/sigma,2));
}

/**
 * 
 */
float pnorm(float q, float mu, float sigma){
  return integrate(dnorm, -100, q, 850, mu, sigma); 
}

/**
 * 
 */
float qnorm(float pvalue, float* z_values, float* p_values, 
            float linspacing, int length_out){
  linfind_t results = linfind_approx(p_values, linspacing, length_out, pvalue); 
  //printf("find:%g, i:%d, count:%d, h:%g\n", results.a0, results.i, results.count, linspacing);
  //printf("find:%g, j:%d, count:%d  h:%g\n", results.a1, results.j, results.count, linspacing);  
  float x0 = p_values[results.i]; float x1 = p_values[results.j];
  float y0 = z_values[results.i]; float y1 = z_values[results.j];
  float q = y0 + (pvalue - x0)*(y1-y0)/(x1-x0);
  return q;
}

void setup() {
  // put your setup code here, to run once:    
  Serial.begin(9600);
  /**
   * Density and Probability function Demo!!!
   */
  int length_out = 50;
  float z_min = -3.00;  float z_max =  3.00;  
  float* z_values = linspace(z_min, z_max, length_out);
  float* p_values = (float*) malloc(length_out * sizeof(float));
  // Compute pvalues from z-values
  for(int i = 0; i < length_out; i++){
      p_values[i] = pnorm(z_values[i], 0.0, 1.0);
  }
  // Print p-values and z-values 
  Serial.println("z-values, p_values");
  Serial.println("------------------");
  for(int i=0; i<length_out; i++){    
    //sprintf(buff, "%f, %f\n", z_values[i], p_values[i]);
    Serial.print(z_values[i], 5);
    Serial.print(", ");
    Serial.print(p_values[i], 5);
    Serial.println();
  }
  //
  // Quantile function using linear interpolation
  //
  float findval = 0.947;
  float h = get_linspacing(-1.00, 1.00, length_out);
  float q = qnorm(findval, z_values, p_values, h, length_out);
  Serial.println();
  Serial.println("Inverse Cumulative Normal Distribution Demo");
  Serial.println("------------------");
  Serial.print("Find Quantile given p-value: "); Serial.print(findval);
  Serial.println();
  Serial.print("q = ");  Serial.println(q, 5);
  //
  // Random Normal Distribution Gen Number Demo 
  //
  Serial.println();
  Serial.println("Random Gaussian Number");
  Serial.println("------------------");
  float rn = 0.0;
  for(int i = 0; i < 100; i++) {
    float r = (float)random(0, 10)/10;
    rn = qnorm(r, z_values, p_values, h, length_out);
    Serial.print(rn); Serial.print(", ");
  }
  //
}

void loop() {
  // put your main code here, to run repeatedly:
}

You need to create a support.h file that only contains the line typedef struct linfind_t linfind_t; in order to download to your atmega.