title: Algorithm for computing the plausibility of a function / Monte Carlo Method
body: <p>I am writing a program that attempts to duplicate the algorithm discussed at the beginning of this article,</p>
  
  <p><a href="http://www-stat.stanford.edu/~cgates/PERSI/papers/MCMCRev.pdf">http://www-stat.stanford.edu/~cgates/PERSI/papers/MCMCRev.pdf</a></p>
  
  <p>F is a function from char to char. Assume that Pl(f) is a 'plausibility' measure of that function. The algorithm is:</p>
  
  <p>Starting with a preliminary guess at the function, say f, and a then new function f* --</p>
  
  <ul>
  <li>Compute Pl(f).</li>
  <li>Change to f* by making a random transposition of the values f assigns to two symbols.</li>
  <li>Compute Pl(f*); if this is larger than Pl(f), accept f*.</li>
  <li>If not, flip a Pl(f)/Pl(f*) coin; if it comes up heads, accept f*.</li>
  <li>If the coin toss comes up tails, stay at f.</li>
  </ul>
  
  <p>I am implementing this using the following code. I'm using c# but tried to make it somewhat more simplified for everyone. If there is a better forum for this please let me know. </p>
  
  <pre><code> var current_f = Initial();    // current accepted function f
   var current_Pl_f = InitialPl();  // current plausibility of accepted function f
  
   for (int i = 0; i &lt; 10000; i++)
          {
              var candidate_f = Transpose(current_f); // create a candidate function
  
              var candidate_Pl_f = ComputePl(candidate_f);  // compute its plausibility
  
              if (candidate_Pl_f &gt; current_Pl_f)            // candidate Pl has improved
              {
                  current_f = candidate_f;            // accept the candidate
                  current_Pl_f = candidate_Pl_f; 
              }
              else                                    // otherwise flip a coin
              {
                  int flip = Flip(); 
  
                  if (flip == 1)                      // heads
                  {
                      current_f = candidate_f;        // accept it anyway
                      current_Pl_f = candidate_Pl_f; 
                  }
                  else if (flip == 0)                 // tails
                  {
                      // what to do here ?
                  }
              }
          }
  </code></pre>
  
  <p>My question is basically whether this looks like the optimal approach to implementing that algorithm. IT seems like I may be getting stuck in some local maxima / local minima despite implementing this method. </p>
  
  <p><b>EDIT</b> - Here is the essentially whats behind Transpose() method. I use a dictionary / hash table of type &lt;&lt; char, char >> that the candidate function(s) use to look any given char -> char transformation. So the transpose method simply swaps two values in the dictionary that dictates the behavior of the function. </p>
  
  <pre><code>    private Dictionary&lt;char, char&gt; Transpose(Dictionary&lt;char, char&gt; map, params int[] indices)
      {
          foreach (var index in indices)
          {
              char target_val = map.ElementAt(index).Value;   // get the value at the index
  
              char target_key = map.ElementAt(index).Key;     // get the key at the index
  
              int _rand = _random.Next(map.Count);   // get a random key (char) to swap with
  
              char rand_key = map.ElementAt(_rand).Key;
  
              char source_val = map[rand_key]; // the value that currently is used by the source of the swap
  
              map[target_key] = source_val; // make the swap
  
              map[rand_key] = target_val;
          }
  
          return map; 
      }
  </code></pre>
  
  <p>Keep in mind that a candidate function that uses the underlying dictionary is basically just:</p>
  
  <pre><code>public char GetChar(char in, Dictionary&lt;char, char&gt; theMap)
  {
       return theMap[char]; 
  }
  </code></pre>
  
  <p>And this is the function that computes Pl(f):</p>
  
  <pre><code>    public decimal ComputePl(Func&lt;char, char&gt; candidate, string encrypted, decimal[][] _matrix)
      {
          decimal product = default(decimal);
  
          for (int i = 0; i &lt; encrypted.Length; i++)
          {
              int j = i + 1;
  
              if (j &gt;= encrypted.Length)
              {
                  break;
              }
  
              char a = candidate(encrypted[i]);
              char b = candidate(encrypted[j]);
  
              int _a = GetIndex(_alphabet, a); // _alphabet is just a string/char[] of all avl chars 
              int _b = GetIndex(_alphabet, b);
  
              decimal _freq = _matrix[_a][_b]; 
  
              if (product == default(decimal))
              {
                  product = _freq;
              }
              else
              {
                  product = product * _freq;
              }
          }
  
          return product;
      }
  </code></pre>