Loop with indexed function not defined for some index values [duplicate]

Question

If anyone cares, this is supposed to be a very simple implementation of a Metropolis algorithm for a uniform distribution. I want to plot the variance of the sampled values as a function of an acceptance ratio, a. The loop in question:

    a = 0.1;
    u[0] = RandomReal[];
    While[a < 2.1,
     k = 1;
     While[k < 11,
      u[k] = RandomReal[];
      ut = u[k] - u[k - 1];
      r = Exp[-a*ut];
      If[r > RandomReal[], k = k + 1, k = k]];
     ui = Table[u[k], {k, 10}];
     m[a] = Sum[u[k], {k, 1, 10}]/10;
     v[a] = Variance[ui];
     a = a + 0.1]

This is working for the most part, however v[a] is undefined for a=0.7, 0.9, and a>1.5. The weird part is that if I check v[0.7] manually it returns a value, but if I make a table with v[a],{a,0.1,2,0.1} it shows as undefined at a=0.7. I want to know why I'm having trouble for these a values.

Answer 1

It is a really bad idea to use floating point numbers for indexing.

In order to see why, have a look at the output of following code after having run your code:

DownValues[v][[16]] // FullForm

HoldPattern[v[1.6000000000000003]] :> 0.10194258571474865

As you can see, it is not v[1.6] which is defined but v[1.6000000000000003]. The deviation of 1.6000000000000003 from 1.6 stems from the fact that you executed a = a + 0.1 repeatedly. Since floating point numbers are stored in binary format, the number 0.1 cannot be expressed exactly, so the tiny error that has to be made may accumulate over time. If the relative error between a number x and each number y for which v has a downvalue is too large, Mathematica won't recognize v[x] as already defined and thus leaves it unevaluated.

So a more robust way is to use integers for indexing:

ClearAll[v];
a = 1;
u[0] = RandomReal[];
While[a < 21, k = 1;
 While[k < 11, u[k] = RandomReal[];
  ut = u[k] - u[k - 1];
  r = Exp[-0.1 a*ut];
  If[r > RandomReal[], k = k + 1, k = k]];
 ui = Table[u[k], {k, 10}];
 m[a] = Sum[u[k], {k, 1, 10}]/10;
 v[a] = Variance[ui];
 a = a + 1
 ]