Project Euler 185

Number MindをMCMCで解いてみた。
方針は与えられた数sに対して、評価関数c(s)を真の正解の数との二乗誤差とおく。今状態sから状態tへの遷移確率を

とおいて遷移していく、これを繰り返していくと

に従ってsをサンプリングすることができ、c(s)が低いものほど頻繁に出やすくなるため、そのうち求める解であるc(s)=0となるsが求められる。

import java.util.*;

public class P185 {
  static final String digits[] = { "5616185650518293", "3847439647293047",
      "5855462940810587", "9742855507068353", "4296849643607543",
      "3174248439465858", "4513559094146117", "7890971548908067",
      "8157356344118483", "2615250744386899", "8690095851526254",
      "6375711915077050", "6913859173121360", "6442889055042768",
      "2321386104303845", "2326509471271448", "5251583379644322",
      "1748270476758276", "4895722652190306", "3041631117224635",
      "1841236454324589", "2659862637316867", };
  static final int correct[] = { 2, 1, 3, 3, 3, 1, 2, 3, 1, 2, 3, 1, 1, 2, 0,
      2, 2, 3, 1, 3, 3, 2 };
  static int cost(int seq[]){
    int cs[] = new int[correct.length];
    for(int i = 0 ; i < correct.length ; i++){
      for(int j = 0 ; j < seq.length ; j++){
        if(seq[j] == digits[i].charAt(j) - '0'){
          cs[i]++;
        }
      }
    }
    int ret = 0;
    for(int i = 0 ; i < cs.length ; i++)
      ret += (cs[i] - correct[i]) * (cs[i] - correct[i]);
    return ret;
  }
  public static void main(String[] args) {
    Random rand = new Random();
    int seq[] = new int[16];
    for(int i = 0 ; i < 16 ; i++)
      seq[i] = rand.nextInt(10);
    int c = cost(seq);
    while(c > 0){
      int ni = rand.nextInt(seq.length);
      int nv = rand.nextInt(10);
      int cv = seq[ni];
      seq[ni] = nv;
      int nc = cost(seq);
      double accept = Math.min(1.0 , Math.exp(- 2.0 * (nc - c)));
      if(rand.nextDouble() < accept){
        c = nc;
      }else{
        seq[ni] = cv;
      }
    }
    System.out.println(Arrays.toString(seq));
  }
}