2

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.

Ryland Forsythe
  • 21
  • 1

1 Answers1

3

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
 ]
Henrik Schumacher
  • 106,770
  • 7
  • 179
  • 309
  • This was the issue, and an easy fix. Thank you – Ryland Forsythe Mar 14 '19 at 19:01
  • You're welcome. Btw.: It may also be a good idea performancewise to define v, u, m as a lists of fixed length (if this length is known before the loop start) and to write into the list with v[[a]] = ... etc. – Henrik Schumacher Mar 14 '19 at 19:03
  • 4
    It is a really bad idea to use floating point numbers for indexing Yes. Actually in Ada (a very well designed and safe language), this will generate a compile error. May be the new Mathematica Lint tool that I saw few days ago in some post can catch such user error, I do not know. – Nasser Mar 14 '19 at 21:20