Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Factor of 2? #6

Open
JonHowMIT opened this issue Nov 5, 2019 · 2 comments
Open

Factor of 2? #6

JonHowMIT opened this issue Nov 5, 2019 · 2 comments

Comments

@JonHowMIT
Copy link

Hi - is is possible that this equation should be modified as noted below?

// return mean + variance * cos( 2 * M_PI * R1) * sqrt(-log(R2));
return mean + variance * cos( 2 * M_PI * R1) * sqrt(-2*log(R2));

This appears to be the Box-Muller Transformation, and thus I think there needs to be a 2 inside the sqrt.
thx

@drf5n
Copy link

drf5n commented Mar 26, 2023

Yes. The Box-Muller transformation does mean that:

return mean + variance * cos( 2 * M_PI * R1) * sqrt(-log(R2));

should read

         return mean + variance * cos( 2 * M_PI * R1) * sqrt(-2*log(R2));

This example of summing samples of Gaussian(0,1) produces approximately mean(0), var=0.5:

// https://wokwi.com/projects/360282019322936321
// for https://github.com/ivanseidel/Gaussian/issues/6

#include "Gaussian.h"


const int sample_size = 10000;
const float var_set = 1.0;
const float mean_set = 0;

void setup() {
  Serial.begin(9600);
  Serial.println("\n\nStarting Random Gaussian Distribution...");
  Serial.print("Gaussian(");
  Serial.print(mean_set);
  Serial.print(",");
  Serial.print(var_set);
  Serial.println(")");
  delay(20);
  randomSeed(20230326);
}


void loop() {
  Gaussian g1 = Gaussian(mean_set, var_set);
  float sum = 0, sumsq = 0, mean = 0, variance = 0;

  // Zero the array
  for (int i = 0; i < sample_size; i++) {
    double x = g1.random();
    sum += x;
    sumsq += x*x;
  }
  mean = sum/sample_size;
  variance = (sumsq-sample_size*mean*mean)/(sample_size-1);

  Serial.print("n:");
  Serial.print(sample_size);
  
  Serial.print(" sum:");
  Serial.print(sum);
  
  Serial.print(" ssq:");
  Serial.print(sumsq);
  
  Serial.print(" mean:");
  Serial.print(mean);
  
  Serial.print(" var:");
  Serial.print(variance);
  Serial.println();
}

produces this output:



Starting Random Gaussian Distribution...
Gaussian(0.00,1.00)
n:10000 sum:-21.59 ssq:5016.92 mean:-0.00 var:0.50
n:10000 sum:-67.49 ssq:4974.67 mean:-0.01 var:0.50
n:10000 sum:54.96 ssq:5090.27 mean:0.01 var:0.51
n:10000 sum:55.73 ssq:4858.84 mean:0.01 var:0.49
n:10000 sum:68.58 ssq:5038.47 mean:0.01 var:0.50
n:10000 sum:-31.55 ssq:4904.44 mean:-0.00 var:0.49
n:10000 sum:56.46 ssq:5057.28 mean:0.01 var:0.51
n:10000 sum:-19.36 ssq:4876.92 mean:-0.00 var:0.49
n:10000 sum:32.28 ssq:4994.59 mean:0.00 var:0.50

See this for a simulation
https://wokwi.com/projects/360282019322936321

image

For Gaussian(0,10) the results are mean=0, var = 50:



Starting Random Gaussian Distribution...
Gaussian(0.00,10.00)
n:10000 sum:-215.91 ssq:501691.40 mean:-0.02 var:50.17
n:10000 sum:-674.85 ssq:497466.50 mean:-0.07 var:49.75
n:10000 sum:549.60 ssq:509027.06 mean:0.05 var:50.90
n:10000 sum:557.26 ssq:485883.06 mean:0.06 var:48.59
n:10000 sum:685.77 ssq:503847.90 mean:0.07 var:50.39
n:10000 sum:-315.49 ssq:490444.37 mean:-0.03 var:49.05
n:10000 sum:564.61 ssq:505727.43 mean:0.06 var:50.57

The input spread parameter appears to be sigma rather than variance.

Gaussian( mu, sigma) versus Gaussian(mu, sigma^2)

drf5n added a commit to drf5n/Gaussian that referenced this issue Mar 26, 2023
@drf5n drf5n mentioned this issue Mar 26, 2023
@drf5n
Copy link

drf5n commented Mar 26, 2023

Looking at it further, the line

https://github.com/ivanseidel/Gaussian/blob/master/Gaussian.h#L130

should be:

return mean + sqrt(variance) * cos( 2 * M_PI * R1) * sqrt(-2*log(R2));

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants