Please weigh in on photon noise

OK, I wrote simpler code from scratch and this time "random(0.000001,1.0)" did not produce an error with Processing's fast floating point math and functions, so when it didn't work completely before, I was probably doing something else, like perhaps adding 0.5 instead of 1.0 before the integer cast, which would output some integer zeros, which makes sense, because I usually got those extra errors (means that suggested a much higher p value than what I input) when the distribution was over a smaller range, with a higher input for p. Adding 0.5 before going integer is normally ideal, but not when you don't want zeros, so in order to have no zeros, you must accept the intervals being slightly larger than they should be. With a p of 0.001 (where the average interval is 1000 units), that should be insignificant, but when p is 0.5, then the intervals are 0.5 units too large on average when most intervals are less than 3 units. Ideally, we'd work with the floating point intervals, but they are more work to visualize with discrete graphical elements.

I wouldn't worry about zero waits in this context as long as the random number generator does its job competently and it is not reset in between calls. After all time is quantized by definition in your approach, so one could argue that when a zero is generated it is simply an artifact of time not being continuous/granular enough - the probability of getting a collision becomes infinitesimally small with an infinitesimally small time period. This is also an unspoken assumption in my approach and the earlier simulation: that Matlab's Mersenne twister, 64-bit double precision random number generator represents continuous time.

In the end imho the litmus test for the script is the resulting distribution. After 20,000 or so arrivals are you getting the expected mean and standard deviation within a fraction of one percent - and the relative Poisson mass function? To see the latter clearly I would suggest histograms with bins that result in a small mean, say around 5 per period. If you choose 5.2 you have my distribution pre-computed for comparison in the link just above.

Jack

That was probably me, agreeing with you that your approach is valid even though technically Binomial rather than straight Poisson. This was my thinking: a number (n) of Bernoulli trials with probability of success (p), otherwise (1-p), is a Bernoulli process, which results in a Binomial distribution; a Binomial distribution as the number of trials approaches infinity and np approaches a finite limit tends to Poisson, with the mean rate equal to np. So according to Wikipedia a rule of thumb says that if n >100 and np<10 you are good to go.

Jack

John Sheehy's approach does generate independent binomial approximations to Poisson.

Jack's code doesn't generate independent, identically distributed, variates:

If we choose a particular slot - say slot 0 - before looking at any of the other slots, the distribution of the histogram count for that slot is B(104000,1/20000) ~ Poisson(5.2) , which, I guess, is what you are trying to achieve.

If we get a count of m in slot 0, that leaves 104000-m counts to put in the other slots.

So if we now choose a second slot - say slot 1 - the distribution of the histogram count for that slot is B(104000-m,1/19999)

After we've looked at 19999 slots, we don't need to look at the 20000th slot to see what the count is there - we can just add up the counts for the 19999 slots that we've already looked at, and take the sum away from 104000.

By fixing the total number of raindrops, the number of raindrops in each bucket becomes dependent on the number of raindrops in the other buckets.

Well, this is where all those degrees with 'mathematics' in them come in handy, so here are more questions. The first Matlab line above generates 104k 64-bit floating point random values between zero and one (then scaled to match the desired experimental range). The interval between zero and one can represent the normalized continuous-time duration of the experiment for a single detector; or the normalized continuous-width of the sensor during a single exposure. Therefore for our purposes the generated random values are specific points in continuous time or space. Aren't these 104k values considered independent, identically distributed variates, assuming a well behaved random number generator?

The second line above simply divides the time/space into a discrete number of contiguous intervals and counts the number of events that happen to fall within each.

Jack

Fixing the total number of arrivals is indeed a limitation of my approach (shades of the earlier Expected Value vs Mean discussion), which however I believe becomes irrelevant with a large enough number of intervals (the dual of JohnS's n and np rule of thumb). What if I told you that I started with a known constant rate of 5.200 arrivals and that I wanted to check 20k such intervals (say a uniformly illuminated 200x100 pixel patch captured in the raw data)?

If the arrivals per interval are independent, and Poisson distributed (identically or not),the total number of arrivals over all intervals is also Poisson distributed, with parameter equal to the sum of the parameters for the individual intervals.

Other kinds of proof are available. For example, along the lines of my earlier description of an arrival process with an arrival rate, and a probability of acceptance.

In the particular case where there are 20000 intervals with distribution of arrivals Poisson(5.2), the total number of arrivals has distribution Poisson(5.2 * 20000 = 104000), or approximately Normally distributed, with mean 104000, and standard deviation sqrt(104000) ~= 332.5 .

The results that you get are at least similar to Poisson, and it intuitively seems likely that you would see that, when you have a lot of raindrops. I have trouble reasoning exactly about this. JACS might do better. More familiar problems are how to divide the raindrops into the buckets. But it's not "self evident" to me that your code makes sense.

I know, that's why these exchanges are always an interesting learning opportunity. I think what makes it not evident at first is that in our universe we always take it as a given that time is monotonic. Looking at the problem in the space domain takes away the monotonicity and makes it easier to understand why it works.

Pretend that light in the shape of a homogeneous slit hits a small portion of a wall in an otherwise dark room, creating a line of normalized length 1 while it is on. Let's say that we allow the source to produce 104,000 photons and then turn it off. These photons arrive at the wall randomly distributed along the length of the line, that is between 0 and 1. Generating 104k homogeneously distributed random numbers between 0 and 1 is the job of the first line of Matlab code. These represent an actual continuous physical position along the line. It does not matter in what order they arrive, what matters is that they are randomly and homogeneously distributed along it.

Next we discover that there are pixels where the line is on the wall, 20k contiguous pixels between normalized position 0 and 1. Each pixel counts the number of incoming photons arriving on its area perfectly. This is the job of the second line of Matlab code. With the provisos mentioned earlier, the distribution of photons counted by the 20k pixels is to an excellent approximation Poisson because the arrivals along the line are many, random, independent and I actually came up with 104k photons because I wanted a mean rate of 5.200 arrivals per interval (er pixel).

It works the same way if we consider a single pixel over 20,000 intervals of time, but monotonicity makes it less intuitive. Would it help the understanding if after generating the 104k random numbers between 0 and 1 I had sorted them? It makes no difference to the end result. Is sorting necessary in quantum mechanics or can we assume that time is not monotonic? Curious.

Jack

You have an ideal single-photon source? 😀

That's not how light sources in photography work.

But, suppose you do have some weird supercooled quantum dot - or whatever - lamp.

So your photodetectors have 100% quantum efficiency, right? And no photoelectrons get lost before readout, and there's no dark current. And no losses due to finite apertures of the emitter or detector.

Eric Fossum has just decided to retire, and spend the rest of his days drinking cocktails on a Pacific beach. 😀😀

What you're describing isn't photography - you're describing an emitter connected to a detector by an optical waveguide. Perhaps part of a quantum computer.

Earlier in this thread, we talked about quantum efficiency, and how to combine that with a Poisson arrival process of photons, giving a special case of a compound Poisson arrival process,because that's modelling something physically realistic.

Now, you want a deterministic arrival process, because your code implies that, rather than because of physical reasoning.

If your magic lamp deterministically ouputs exactly one photon per picosecond into a lossless optical waveguide, but the photodetector has some quantum efficiency, 0 < p < 1, John Sheehy's code models exactly that.

If you use your magic 1 photon/ps lamp to illuminate a real scene, and capture photons from the scene with a camera with finite aperture(s), and typical QE, the ratio of detected photons to emitted photons will be << 1. Also, the arrival times of the detected photons will have varying delays relative to the emitted photons.

Your model departs from physical plausibility at the point at which you say you know how many photons are emitted.

Your code does indeed produce approximately Poisson-distributed variates. But you seem to be attempting to claim some physical significance for your code from this. Which is circular reasoning.

If you just want to generate some Poisson-distributed variates in matlab/octave, just use poissrnd(...). That will do a much better job.

You need to start with a physically plausible model, and then show that it results in a Poisson distribution.

What is “homogeneously distributed” random numbers?

I expect he means "uniformly distributed".

My last few posts (and I thought yours and JohnS's) were not about the detection of light, we discussed that interaction earlier; they were about showing via Monte Carlo simulation that random, uniformly distributed cars/photons arriving at a certain mean rate, at a given point in space/time, result in a Poisson distribution (measured by a perfect detector for simplicity). I was pleasantly surprised that it was so easy to do in two intuitive, if not self-evident, lines of code. In fact originally there were three lines because I sorted the random generated values before binning them - but then I realized that sorting was not necessary.

Not so much emitted as arriving: if we have the rate and the number of intervals we know how many cars/photons are expected, right? I believe knowing the exact value ex-post or ex-ante are two sides of the same coin, it was the earlier point of discussion between Bob and JACS. Imho there is no getting around a certain amount of uncertainty either way, it's just that the uncertainty affects different variables in each case.

In photography there is a finite Exposure time so we are constrained by that and cannot wait till the 104,000th photon has arrived before stopping the counting. On the other hand for a given Exposure time we can look at how many random, uniformly distributed photons have arrived on the area covered by our putative perfect pixels. From there we can compute an expected value, the more the samples the less the uncertainty, which is generally true.

Rather than a claim it's curiosity. Is there something in QM that says that time needs to increase monotonically?

The idea that, say, a huge, nearby, Lambertian blackbody radiator like the sun only produces one photon at a time in our direction boggles the mind. Then again, QM often does.

Jack

"Increase" and "monotonically" say the same thing here but increase with respect to what?

The Poison process assumptions do not really require that. First, the relevant assumption is about a probability; and then you need to take a small time interval. That assumption then is that the probability of two or more particles over that interval decreases faster than the length of the interval.

I guess the naive question is whether we can assume that all photons could be emitted randomly timelessly - and that it is observers like us who live in a universe subject to time as we know it that are forced to impose a monotonically increasing (as opposed to decreasing;-) structure on them.

Makes sense, thanks JACS.
Jack

This belongs to relativity. There is relativistic QM, of which I know nothing. The appropriate model here would be just relativity though. Then photons are just “particles” moving with speed of light. More precisely, they are null geodesics of the background Minkowski, or more generally, Lorentzian metric. Now, there is no “god given” time in relativity. Not even “god given” time direction like past and future - you can reverse a fixed time orientation, and “nothing much changes” in the theory. Locally, you can choose many time functions, indeed, but their gradients must be in one of the two timelike cones; one of which we can call future oriented, the other one: past oriented. Once we made that choice at the source, and emit a photon in the so-defined future, that photon “remembers” its orientation. Say it arrives at us much later, like from the Big Bang to today. We are free to choose time orientation now as we wish but only one of those two choices would be consistent with the future orientation of that photon which it carries with itself.

This is getting a bit obscure and OT, but I wanted to tie up a couple of loose ends about generating geometrically distributed variates.

I may have been a bit dim in replying earlier.

Geometrically distributed variates are distributed over the non-negative integers. Whis is quite a big range. Much larger than "int", "long", or even "double".
So there's a bit of a problem about what to do with extreme values, which become very likely when p is close to the (positive) lower end of the "double" range. We have an integer-valued distribution, but variates may not fit in a "long", when p is small. We can just clamp them to MAX_LONG, but that will inevitably distort the distribution, particularly for sufficiently small values of p. But very small values of p imply variates that overflow the "double" range, never mind "long".

I've included some C++ code for generating geometrically distributed variates, which should translate easily to Java. It might seem a bit long, but it's mostly comments. Should compile with recent g++, giving an executable file. There's two functions for generating geometrically distributed variates: "geometric_b()", and "geometric_e()". The second one is a bit faster sometimes.

There's also a couple of sequence generators in there. The important one is the "golden" sequence generator, which is useful for 1D sampling problems generally. Starting with the interval [0,1), successive outputs from the sequence generator divide the largest previous subinterval by the golden ratio.

#include <cstdlib>
#include <cstdio>
#include <climits>
#include <cmath>
#include <random>
#include <iostream>
#include <string>

using namespace std;

mt19937_64 RNG; // 64-bit Mersenne twister Random Number Generator
uniform_real_distribution<double> u01(0.0,1.0);


// Sequence generator for testing distribution generators:
// Outputs the largest doubles < 1, in decreasing order.
// Reaches 0.0 after 2^53 - 1 steps, then cycles back to (double)(2^53 - 1)/(2^53):

double high_sequence()
{
  const double one_over_two_to_the_53 = 1.0 / ( (double)(1 << 30) * (double)(1 << 23) );
  const long mask = (1l << 53) - 1;

  static long state = mask;
  double ret = (double) state * one_over_two_to_the_53;
  // printf( "high: 0x%.14lx %.17f\n", state, ret );
  state = (state - 1) & mask;
  return ret;
}

// Sequence generator for testing distribution generators:
// Outputs a quasirandom sequence that uses properties of the golden ratio
// to explore [0,1) as evenly and quickly as possible (with minimum "discrepancy").

double golden_sequence()
{
  const double one_over_two_to_the_53 = 1.0 / ( (double)(1 << 30) * (double)(1 << 23) );
  const long mask = (1l << 53) - 1;

  // Odd number nearest to 2^53 * (sqrt(5)-1)/2 (computed to arbitrary precision):
  const long golden = 0b10011110001101110111100110111001011111110100101001111;

  static long state = 0;
  double ret = (double) state * one_over_two_to_the_53;
  // printf( "golden: %.17f\n", ret );
  state = (state + golden) & mask;
  return ret;
}

double U;   // Random variate in [0,1)


// Relatively straightforward generator for (0-based) geometrically distributed variates.
// log-based, but with an initial cull, which causes execution time to fall as
// p increases.

long geometric0_b( double p )
{
  // double U = u01(RNG);   // u01 guarantees 0 <= U < 1
  double U = ::U;

  // if U < p, the waiting time is zero, and we take a quick exit.
  //
  // This test also culls cases p >= 1, and p is a NaN.

  if( ! (U >= p) )  // IEEE floating-point comparison.
    return 0;

  // We have:
  //  p is not NaN
  //  U is not NaN
  //  p <= U < 1

  // p == 0.0 implies an infinite wait.
  //
  // For sufficiently small p > 0, the wait is greater than LONG_MAX, regardless of U.
  // (We have already culled the case 0  <= U  < p).
  // The smallest possible value of U which is > 0 is 2^-53
  //
  // A test is needed here to cull p <= 0, so we might as well cull some more cases
  // with the test while we're here.

  double overflow_threshold = 1.0 / ( (double) (1l << 53) * (double) LONG_MAX );  // 2^-116

  if( p < overflow_threshold ) 
    return LONG_MAX;

  // We have: overflow_threshold <= p <= U < 1.0

  // The following expression may have a value greater than LONG_MAX,
  // if U is sufficiently large, and p is sufficiently small, but cannot overflow
  // the range of normalised double (because of earlier culling).

  double z = log1p(-U) / log1p(-p);

  // In C/C++ the result of converting a floating point value to an integer type is undefined
  // when the floating point value is greater than the largest value of the integer type.
  // Some machines will give LONG_MAX here; some will give garbage.
  // This test is not needed in Java (Java guarantees LONG_MAX).

  // LONG_MAX is not exactly representable as a double on systems where long is > 54 bits.
  // Such as Java, and LP64 systems.
  // On these systems, LONG_MAX rounds to a double value which is 1.0 greater than the long value LONG_MAX.

  if( ! (z < (double) LONG_MAX) )   // IEEE floating-point comparison.  z should not be NaN, but...
    return LONG_MAX;

  return (long) z;  // Truncates towards zero.
}

// Generator for (0-based) geometrically distributed variates.
// Uses either log-based, or Cumulative Distribution Function based probability estimation,
// mainly depending on the value of p.
// Speed increases with p.   Significantly faster than simpler generator above when p >= 1/2,
// but p < about 0.9

long geometric0_e( double p )
{
  // double U = u01(RNG);
  double U = ::U;

  // Assume that U >= 0 and U < 1 are guaranteed.

  // if U < p, the waiting time is zero, and we take a quick exit.
  //
  // This test also culls cases p >= 1, and p is a NAN.

  if( ! (U >= p) )  // IEEE floating-point comparison.
    return 0;

  // We have:
  //  p is not NaN
  //  U is not NaN
  //  p <= U < 1
  //  0 <= U

  // Choose between logarithmic and Cumulative Distribution methods,
  // depending on the values of p and U:
  //
  // If p < 1/2, always use the logarithmic method. This avoids the potential of huge numbers
  // of iterations of the CDF method, with associated accumulation of rounding errors.
  //
  // If p >= 1/2, the CDF method is *usually* more accurate, and usually faster.
  // However, for values of p slightly above 1/2, and values of U quite close to 1,
  // the CDF method can require tens of iterations, losing accuracy and speed.
  //
  // Some of these problems could be resolved by creating a class that generates
  // geometric variates for fixed p, with a constructor that chooses between different
  // generation methods based on the value of p.

  double force_log_threshold = 255.0/256.0;

  if( p < (1.0/2.0) || U >= force_log_threshold )
  {
    // p == 0.0 implies an infinite wait.
    //
    // For sufficiently small p > 0, the wait is greater than LONG_MAX, regardless of U.
    // (We have already culled the case 0  <= U  < p).
    // The smallest possible value of U which is > 0 is 2^-53
    //
    // A test is needed here to cull p <= 0, so we might as well cull some more cases
    // with the test while we're here.

    double overflow_threshold = 1.0 / ( (double) (1l << 53) * (double) LONG_MAX );  // 2^-116

    if( p < overflow_threshold ) 
      return LONG_MAX;

    // We have: overflow_threshold <= p <= U < 1.0

    // The following expression may have a value greater than LONG_MAX,
    // if U is sufficiently large, and p is sufficiently small, but cannot overflow
    // the range of normalised double (because of earlier culling).

    double z = log1p(-U) / log1p(-p);

    // In C/C++ the result of converting a floating point value to an integer type is undefined
    // when the floating point value is greater than the largest value of the integer type.
    // Some machines will give LONG_MAX here; some will give garbage.
    // This test is not needed in Java (Java guarantees LONG_MAX).

    // LONG_MAX is not exactly representable as a double on systems where long is > 54 bits.
    // Such as Java, and LP64 systems.
    // On these systems, LONG_MAX rounds to a double value which is 1.0 greater than the long value LONG_MAX.

    if( ! (z < (double) LONG_MAX) )
      return LONG_MAX;

    return (long) z;    // Truncates towards zero.
  }

  // p >= 1/2, and we usually get more accurate probabilities
  // by computing the cumulative probabilities directly.

  // We have: 1/2 <= p <= U < force_log_threshold
  // Which implies U-p is computed exactly. (At least in the first loop iteration).

  double q = 1.0 - p;   // q is exactly representable - no rounding error in computing q.
  long n = 0;

  do
  {
    n = n + 1;
    U = U - p;  // Subtracting current probability from U gives less rounding error than accumulating the CDF.
    p = p * q;  // p <= 2^-(n+1).  Maximum error accumulation in p <= 0.5 ulp / iteration.
  }
  while( U >= p );

  return n;
}

int main(int argc, char ** argv)
{
  double p = argc > 1 ? atof(argv[1]) : 1.0/2.0;

  printf( "p  =  %.17f\n", p );

  const int n = 500000000;

  double x_max = 0;
  double mean = 0;
  double sx = 0;
  double s_var = 0;

  double max_gap = ( -53 * log(2) / log1p(-p) );

  for( int i = 0; i < n; ++i )
  {
    // U = u01(RNG);
    // U = high_sequence();
    U = golden_sequence();
#if 0
    // printf( "U = %.16f\n", U );
    long b = geometric0_b(p);
    long e = geometric0_e(p);
    // printf( "geom_b =%3ld\n", b );
    // printf( "geom_e =%3ld\n", e );

    if( (b != e) || (b > max_gap) )
    {
      printf( "Ouch!\n" );
      printf( "U = %.17f; i=%d\n", U, i );
      printf( "Unrounded log estimate: %.17f\n", log1p(-U) / log1p(-p) );
      printf( "geom_b =%3ld\n", b );
      printf( "geom_e =%3ld\n", e );
      printf( "max_gap = %3f\n", max_gap );
      return 1;
    }
#else
    long e = geometric0_e(p);
#endif

    if( e > x_max ) x_max = e;

    double emm1 = e - mean; // Welford's method for stable variance.
    sx += e;
    mean = sx / (double) (i+1);
    s_var += emm1 * (e - mean);
  }

  printf( "Max = %9g; Max possible = %9g; max/max_possible = %.5f\n", x_max, max_gap, x_max/max_gap );
  printf( "Mean     %16g; expected %16g\n", mean, (1-p) / p );
  printf( "Variance %16g; expected %16g\n", s_var / n, (1-p) / (p*p) );

  return 0;
}