Modelica.Math.Random.Utilities.impureRandomInteger: Not an uniform distribution for imin=1

[jmoeckel]

The algorithm of this function is:

  r  := impureRandom(id=id);
  y  := integer(r*imax) + integer((1-r)*imin);
  y  := min(imax, max(imin, y));

where y is the output and should be (as stated in the documentation) "A random number with a uniform distribution on the interval [imin,imax]".

This is not true if imin=1. In this case, the part integer((1-r)*imin) will always be 0 as r<=1. Therefore, the first computation of y results in something smaller than imax except if r=1. Given the second equation for y this is also the mandatory condition for y=imax, while the probability for all other integers on the intervall [imin, imax] is much higher (actually I think, that the probability for imin is twice than the one for any integer between imin and imax because of the max(imin, y) part).

I do not think, that this is something to call a "uniform distribution".

[HansOlsson]

That seems like an incorrect implementation and we should instead use something like:
y:=min(imax, integer(r*(imax-imin+1))+imin);
The min is only used in the very rare case that r=1.

Testing the formula with imin=0 and imax=N gives that the value before truncation is almost [0,N+1) and thus after truncation each integer has almost 1/(N+1) probability.

We could also use the normal name:
https://en.wikipedia.org/wiki/Discrete_uniform_distribution

[beutlich]

Test model to play with:

model T2252
  parameter Integer id(fixed=false);
  Integer y;
  initial algorithm
    id := Modelica.Math.Random.Utilities.initializeImpureRandom(123456789);
  equation
    when sample(0, 0.001) then
      y = Modelica.Math.Random.Utilities.impureRandomInteger(id, 1, 3);
    end when;
end T2252;

[HansOlsson]

I updated the test-model to be test that the values are reasonable - and make it possible to generalize:

model T2252
  parameter Integer id(fixed=false);
  Integer y;
  parameter Integer n=3;
  Integer cnt[n];
  parameter Real nrSigma=3;
  Real avg=sum(cnt)/(n);
  Real lowBound=avg-nrSigma*sqrt(avg);
  Real highBound=avg+nrSigma*sqrt(avg);
initial algorithm 
    id := Modelica.Math.Random.Utilities.initializeImpureRandom(123456789);
equation 
    when sample(0, 0.001) then
      y = Modelica.Math.Random.Utilities.impureRandomInteger(id, 1, n);
      cnt=pre(cnt)+{if i==y then 1 else 0 for i in 1:n};
    end when;
    when terminal() then
      for i in 1:n loop
        assert(cnt[i]>lowBound,
          "Number of generated "+String(i)+" is "+String(cnt[i])+" but should be "+String(avg)+"+/-"+String(sqrt(avg)));
        assert(cnt[i]<highBound,
          "Number of generated "+String(i)+" is "+String(cnt[i])+" but should be "+String(avg)+"+/-"+String(sqrt(avg)));
      end for;
    end when;
end T2252;

[beutlich]

@HansOlsson Can you pls manage to open a pull request against master branch for the fix and the test model for ModelicaTest. Thanks a lot.

[beutlich]

Resolved in master. Thanks.