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.
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.
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.
"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.
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;
}
