読者です 読者をやめる 読者になる 読者になる

任意の画像を数式で表現する

WolframAlphaで「Popular Curves」などと検索をかけると、ロゴや食品をはじめ、アニメキャラクターや実在する人物までを媒介変数表示された曲線で表している作品が多く見受けられる。下の画像はアインシュタインのプロットらしい。

f:id:I-was-a-Ki:20170122161926p:plain

このように、任意の画像を数式で表現する方法については、以下のリンクに詳しい記載がある。 aomoriringo.hateblo.jp

しかし、上記のサイトは有償ソフトウェアたるMathematicaを利用しているためハードルが高い。そこで、Mathematicaを購入することなく利用する方法を解説した下記のサイトを参考にし、両者を組み合わせることによって画像の数式表現を実現したログを残しておこうと思う。

hirax.net::「Wolfram CDF PlayerをMathematicaとして使う方法」をRubyでもっと簡単にしてみた

なお、以下の作業はWindows10上で行っており, 全体を通じてブラウザはFirefoxを利用している。他のブラウザを使用していて動かない場合は変更を検討してみると良いかもしれない。

A. Mathematicaを利用できる環境を整える

Mathematica自体は有料だが、Mathematicaにより生成されたCDFファイルを再生する「wolfram CDF player」は無料配布されている(MathematicaとCDF playerの関係についてはこちらを参照願いたい)。今回はこれとRubyのWebサーバーを組み合わせることにより実行環境を構築することとする。

wolfram CDF playerの入手

まずは以下のURLからこれを入手しよう。ダウンロード・インストールにそれなりに時間がかかるので注意。

www.wolfram.com

Ruby実行環境の構築

こちらを参考にRubyのセットアップをしよう。以下のURLからRubyインストーラをダウンロードして実行すれば簡単にインストールすることができる。個人的にはパスも通しておいた方が何かと便利なのではないかと思う。

RubyInstaller for Windows

それでは、早速冒頭で述べた記事に記載のrubyコードをコピペして適当なディレクトリに"math.rb"などどして保存し、コマンドプロンプト

ruby math.rb

と実行してみる。Windows10ではここで警告が出たが適当にOKしておこう。なお、このサーバーを終了させるにはhttp://localhost/shutdownにアクセスすればよい(ブラウザにお気に入り登録しておくと便利かもしれない)。

長い文字列を送信できるようにしておく

以下の作業で、WEBから一定以上の数の文字列を送信した場合に「ERROR WEBrick::HTTPStatus::RequestURITooLarge」というエラーメッセージが生じて送信が失敗する場合がある。その場合、こちらを参考に、C:\Ruby23-x64\lib\ruby\2.3.0\webrickに存在するhttprequest.rbをテキストエディタで開き、MAX_URI_LENGTHを指定している箇所を

MAX_URI_LENGTH = 64 * 1024 * 1024 # :nodoc:

などと変更し上書き保存することでこのエラーに備えておくことができる。

Universal Mathematica Manupulator 3 にアクセス

上のRuby Webサーバーを立ち上げた状態で、ブラウザでUniversal Mathematica Manipulator 3にアクセスする。初回のアクセス時に「umm3.cdfをダウンロードしますか」というようなことを訊かれるので適当なディレクトリにダウンロードしておく。上手く行っていれば暫く待っていると下図のような立体が右側領域に表示されるはずである。

f:id:I-was-a-Ki:20170122151320p:plain

※「クラッシュしました」といったメッセージが時たま表示されるが、ブラウザやあるいはPCごと再起動してみたりを何回か行えば(?)復旧する。この辺のことは全然分からないので詳しい方はぜひ教えてほしい。

使い方としては、まず左側に数式を打ち込んだ状態でTABキーを押して(入力内容がRubyのWebサーバーに送信される)、その後右側領域のevalボタンを押せば出力が表示される、という流れである。例えば左側に「1+1」を入力して上記の作業を行った場合、右側には2が表示される。

B. 実際に画像を数式化してみる

注:以下に示すコードは、これも最初に述べた記事に載っているものを一部改変しただけでほとんどコピペである。
今回は、主に筆者がスピッツファンだという理由によりスピッツのロゴを数式で表現してみようと思う。画像のローカルのパスを渡すことによりインポートすることも少なくともデスクトップ版Mathematicaでは可能なようだが、今回は上手くいかなかった(?)

f:id:I-was-a-Ki:20170122153827j:plain:w300

関数の読み込みと画像の取り込み

最初に指定されているmaxOrderNumが使用する数式の最大次数となる。

(* parameters *)
(* 最大次数 *)
maxOrderNum = 200;
(* 画像URL, もしくはローカルパス *)
imageURL = "https://yt3.ggpht.com/-wbOANd5n5ZY/AAAAAAAAAAI/AAAAAAAAAAA/ELOj6QE_dyc/s900-c-k-no-mo-rj-c0xffffff/photo.jpg";

pointListToLines[pointList_, neighbothoodSize_: 6] :=
 Module[{L = DeleteDuplicates[pointList], NF, lambda,
   lineBag, counter, seenQ, sLB, nearest,
   nearest1, nextPoint, couldReverseQ, d, n, s},
  NF = Nearest[L];
  lambda = Length[L];
  Monitor[
   (*list of segments *)
   lineBag = {};
   counter = 0;
   While[counter < lambda,
    (*new segment*)
    sLB = {RandomChoice[DeleteCases[L, _?seenQ]]};
    seenQ[sLB[[1]]] = True;
    counter++;
    couldReverseQ = True;
    (*complete segment*)
    While[
     (nearest = NF[Last[sLB], {Infinity, neighbothoodSize}];
      nearest1 = 
       SortBy[DeleteCases[nearest, _?seenQ], 
        1. EuclideanDistance[Last[sLB], #] &];
      nearest1 =!= {} || couldReverseQ),
     If[
      nearest1 === {},
      (*extend the other end; penalize sharp edges*)
      sLB = Reverse[sLB]; 
      couldReverseQ = False,
      
      (* prefer straight continuation *)
      nextPoint = If[Length[sLB] <= 3,
        nearest1[[1]],
        d = 1. Normalize[(sLB[[-1]] - sLB[[-2]]) +
            1/2 (sLB[[-2]] - sLB[[-3]])];
        n = {-1, 1} Reverse[d];
        s = Sort[{Sqrt[(d.(# - sLB[[-1]]))^2 +
               (*perpendicular *)2
                (n. (# - sLB[[-1]]))^2], #} & /@ nearest1];
        s[[1, 2]]];
      AppendTo[sLB, nextPoint];
      seenQ[nextPoint] = True;
      counter++]];
    AppendTo[lineBag, sLB]];
   (*return segments sorted by length*)
   Reverse[SortBy[Select[lineBag, Length[#] > 12 &],
     Length]],
   (*monitor progress*)
   Grid[
    {{Text[Style["progress point joining",
        Darker[Green, 0.66]]],
      ProgressIndicator[counter/lambda]},
     {Text[Style["number of segments",
        Darker[Green, 0.66]]],
      Length[lineBag] + 1}},
    Alignment -> Left, Dividers -> Center]]]
 
(* Fourier coefficients of a single curve *)
fourierComponentData[pointList_, nMax_, op_] :=
 Module[{epsilon = 10^-3, myu = 2^14, M = 10000, s, scale, delta, L, 
   nds, sMax, if, xyFunction, X, Y, XFT, YFT, type},
  (* prepare curve *)
  scale = 
   1. Mean[Table[
      Max[fl /@ pointList] - 
       Min[fl /@ pointList], {fl, {First, Last}}]];
  delta = EuclideanDistance[First[pointList], Last[pointList]];
  L = Which[
    op === "Closed",
    type = "Closed";
    If[First[pointList] === Last[pointList], pointList,
     Append[pointList, First[pointList]]],
    
    op === "Open",
    type = "Open";
    pointList,
    
    delta == 0.,
    type = "Closed";
    pointList,
    
    delta/scale < op,
    type = "Closed";
    Append[pointList, First[pointList]],
    
    True,
    type = "Open";
    Join[pointList, Rest[Reverse[pointList]]]];
  (*re-parametrize curve by arclength *)
  xyFunction = BSplineFunction[L, SplineDegree -> 4];
  nds = NDSolve[
    {s'[t] == Sqrt[xyFunction'[t].xyFunction'[t]],
     s[0] == 0}, s, {t, 0, 1}, MaxSteps -> 10^5, PrecisionGoal -> 4];
  (* total curve length *)
  sMax = s[1] /. nds[[1]];
  if = Interpolation[
    Table[{s[rho] /. nds[[1]], rho}, {rho, 0, 1, 1/M}]];
  X[t_Real] := 
   BSplineFunction[L][Max[Min[1, if[(t + Pi)/(2 Pi) sMax]], 0]][[1]];
  Y[t_Real] := 
   BSplineFunction[L][Max[Min[1, if[(t + Pi)/(2 Pi) sMax]], 0]][[2]];
  (* extract Fourier coefficients *)
  {XFT, YFT} = 
   Fourier[Table[#[N@t], {t, -Pi + epsilon, 
        Pi - epsilon, (2 Pi - 2 epsilon)/myu}]] & /@ {X, Y};
  {type, 2 Pi/
     Sqrt[myu]*((Transpose[
         Table[{Re[#], Im[#]} &[Exp[I k Pi] #[[k + 1]]], {k, 0, 
           nMax}]] & /@ {XFT, YFT}))}]
Options[fourierComponents] =
  {"MaxOrder" -> maxOrderNum, "OpenClose" -> 0.025};
 
fourierComponents[pointLists_, OptionsPattern[]] :=
 Monitor[
   Table[fourierComponentData[
     pointLists[[k]],
     If[Head[#] === List, #[[k]], #] &[OptionValue["MaxOrder"]],
     If[Head[#] === List, #[[k]], #] &[OptionValue["OpenClose"]]
     ], {k, Length[pointLists]}],
   Grid[
    {{Text[
       Style[
        "progress calculating Fourier coefficients",
        Darker[Green, 0.66]]],
      ProgressIndicator[k/Length[pointLists]]}},
    Alignment -> Left, Dividers -> Center]] /; Depth[pointLists] === 4
 
makeFourierSeries[
  {"Closed" | "Open", {{cax_, sax_}, {cay_, say_}}},
  t_, n_] :=
 {Sum[If[k == 0, 1/2, 1] cax[[k + 1]] Cos[k t] + 
    sax[[k + 1]] Sin[k t], {k, 0, Min[n, Length[cax]]}],
  Sum[If[k == 0, 1/2, 1] cay[[k + 1]] Cos[k t] + 
    say[[k + 1]] Sin[k t], {k, 0, Min[n, Length[cay]]}]}
 
(* 画像読み込み *)
img = Import[imageURL];

これをAで準備したUniversal Mathematica Manipulator 3の左側領域にコピペして実行する(NULLと表示されるはずである)。

エッジを抽出し線分の数を表示

(* エッジ抽出 *)
edgeImage = Thinning[EdgeDetect[ColorConvert[
     ImagePad[Image[Map[Most, ImageData[img], {2}]], 20, White], 
     "Grayscale"]]];
edgePoints = {#2, -#1} & @@@ Position[ImageData[edgeImage], 1, {2}];
SeedRandom[2];
hLines = pointListToLines[edgePoints, 16];

(* 線分の数を表示 *)
Length[hLines]

これを実行し暫く待つと右側に線分の数が表示される。今回使用したスピッツのロゴの場合は3であった。

エッジ抽出の結果を描画してみる

Graphics[{ColorData["DarkRainbow"][RandomReal[]],Line[#]}&/@hLines]

実行結果は下図のようになった。

f:id:I-was-a-Ki:20170122155433p:plain

フーリエ変換を実行&数式を表示

fCs = fourierComponents[hLines];

curves = makeFourierSeries[#, t, 200] & /@ fCs;
curves2 = Rationalize[curves, 0.002];
Style[curves2, 6] // TraditionalForm

少し時間はかかるが、これでめでたく数式を得ることに成功する。 得られた数式はコピペでテキストファイルなどに保存することが可能である。

f:id:I-was-a-Ki:20170122155751p:plain

将棋:一致率の統計解析

遅ればせながらあけましておめでとうございます。
さて、今回は前回の続きで、プロの将棋の一致率計算結果をもとにいろいろと計算してみるとする。ここでは、

  • 一致率の高い方が勝つ確率

  • 手数と一致率の相関関係

をみていくことにしよう。データは2chkifuの00001-01000(番号が飛び飛びになっているため実際は955局しかないようだ)を初手から一手一秒の技巧で解析したもの。詳細は前の記事を参照されたい。

一致率の高い方が勝つ確率

解析結果(ここではresults.txt)にsetwdし、読み込む。

data <- read.table("result.txt")
gote = 0 #一致率が後手>先手かつ後手勝ちの結果を格納する
sente = 0 #一致率が先手>後手かつ先手勝ちの結果を格納する
for (i in 1:length(data[[1]])){
    if(data[[7]][i] %% 2 == 0 && data[[6]][i] > data[[5]][i]){ 
        gote <- gote + 1
    }
    if(data[[7]][i] %% 2 == 1 && data[[6]][i] < data[[5]][i]){
        sente <- sente + 1
    }
}

 

結果は以下のようになる。

> sente
[1] 389
> gote
[1] 292
> (gote + sente) / length(data[[1]])
[1] 0.713089

以上より、一致率の高い方が勝つ確率は約71%と計算される。

手数と一致率の相関関係

ここでは勝った方の一致率を採用して、手数との一致率の相関関係を研鑽してみる。

MR <- rep(0,length(data[[1]]))

for (i in 1:length(data[[1]])){
    if(data[[7]][i] %% 2 == 0){
        MR[i] <- data[[6]][i]
    }
    if(data[[7]][i] %% 2 == 1){
        MR[i] <- data[[5]][i]
    }
}

 

結果は以下のようになる。

> cor(data[[7]],MR)
[1] -0.1929392
> cor.test(MR,data[[7]])

        Pearson's product-moment correlation

data:  MR and data[[7]]
t = -6.0702, df = 953, p-value = 1.842e-09
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.2532767 -0.1311064
sample estimates:
       cor
-0.1929392

以上より相関係数は-0.19と計算され、手数と一致率の間に相関関係はないことが判明した。今回は思考停止して初手から一致率を計算しているので、仕掛け以降の手数を対象にしたら異なる結果が得られる可能性は大いにあろう。

将棋:一致率一括解析ログ

今回の三浦九段にまつわる一連の騒動を受け、(それなりに out of dateではあるが)プロ棋士の将棋の一致率解析を自分でやってみようと思い立った方のためのメモ。因みに将棋のエンジン自体はすでにインストールされているものとし、ここでは技巧を用いる。

棋譜のダウンロード

今回は以下のリンクから2chkifuのデータをダウンロードした。きちんと確認してはいないが1950年代のものから2000年代初頭のものまで合計約5万局、一部女流のものも含まれているようだ。全て展開すると5~10分ほどかかるので注意。
https://code.google.com/archive/p/zipkifubrowser/downloads

②KI2ファイルをKIFファイルに変換

のちのち使用するスクリプトがKIFファイルの解析用なのでこの作業を行う。以下のリンクからKKC.exeをインストールして実行、変換元と変換先のディレクトリを指定すれば一括変換が可能。
http://www.geocities.jp/shogi_depot/

③一致率解析をするスクリプトを入手

以下のComputer Shogi File Uploaderから「USI将棋エンジン一致率計算プログラム usi.vbs Ver.1.2」をダウンロード。
http://mucho.girly.jp/cgi/shogiup/upload.html

④得たいデータに応じて出力部分を書き換える

今回は試しに下のようなデータ(「先手名前」「後手名前」「平均NPS」「先手狭義一致率」「後手狭義一致率」「手数」)を得たいとする。

先手:羽生善治 後手:谷川浩司 738 52.4 49.2 126
先手:島  朗 後手:三浦弘行 633 69.1 76.4 110
先手:羽生善治 後手:谷川浩司 748 50.8 52.4 126 64
先手:島  朗 後手:三浦弘行 651 70.9 78.2 110 56
先手:羽生善治 後手:谷川浩司 746 52.4 49.2 126 

先程DLしたusi.vbsを適当なテキストエディタで開き、これを出力するように書き換えてみよう。

技巧の追加

このスクリプトはもともとGPSShogi/GPSfish/Bonanzaのみ対応のようなので、70行目付近の「思考条件を送る(エンジン別)」の箇所に以下を追加して技巧を追加する。

If     Instr(LCase(strEngine), "gikou") > 0 Then
        strEngineName = "gikou"
        strPonder = "setoption name USI_Ponder value false"
        strThread = "setoption name Thread value " & intThread

出力部分の書き換え

  • ヘッダーの削除:110行目付近にある「ヘッダー出力」の部分の行頭に「'」を入れてコメントアウトする。

  • 改行しないようにする:元ファイルの283-286,292-295行目にある改行の部分を同様にコメントアウト

  • 元ファイル297行目からの「結果まとめ出力」を次のように書き換える。

'結果まとめ出力
str1 = strName(0)
str2 = strName(1)
If numNPS > 0 Then
    str3 = CLng(sumNPS / (1000 * numNPS))
Else
    str3 = ""
End If
str4 = FormatNumber(100.0 * numMatch(0) / numTesu(0), 1)
str5 = FormatNumber(100.0 * numMatch(1) / numTesu(1), 1)
str6 = numTesu(0) + numTesu(1)
    
WriteLog objLog, str1 & " " & str2 & " " & str3 & " " & str4 & " " & str5 & " " & str6 & vbCRLF

その下のファイル区切り部分もコメントアウトする。

⑤実行

④をgikou.exeが存在するディレクトリに保存し(”usi2.vbs”とした)、コマンドプロンプトでそのディレクトリに移動し以下を実行する。以下は一例で、C:\Users\ユーザー名\Downloads\2chkifu\experimentにあるkifファイルを、思考時間1秒4スレッドで解析する。

C:\Users\ユーザー名\Downloads\gikou_win_20160606\gikou_win_20160606>cscript usi2.vbs gikou.exe C:\Users\ユーザー名\Downloads\2chkifu\experiment 4 1

 

これにより、同じディレクトリに「解析結果.txt」という、④の冒頭に記したようなファイルが生成される。次回はこれを使って遊んでみるとしよう。

将棋ソフト騒動の統計的考察 - ベイズ流アプローチのログ -

2016年10月12日に発生した三浦九段竜王戦出場停止の件に関して一致率の解析を行った記録。同様の試みは既出だが(下記リンクを参考にさせて頂いた)、統計の勉強を兼ねて有意差検定ではなくベイズ統計を用いて結果を求めてみようと思う。

design.syofuso.com

注釈:中の人について。古典的な統計は一般教養で単位を取得しただけ、ベイズ統計は入門書を一冊読んだだけという素人である。用語が違う、コードが汚い、そもそもやっていることがおかしい、その他諸々問題点を発見された方はご教授願いたい。

概要

  • 使用したデータ
    以下のサイトより、仕掛けから終局までで求めた一致率を借用。如何なるデータを用いるかについて、恣意的な要素が介在しうることの危険性は議論に値するがここでは深入りしないこととする。
    i2chmeijin.blog.fc2.com

  • 解析方法
    三浦九段の6月以前の対局の一致率を第1群、7月以降の対局の一致率を第2群として扱う。初めにそれぞれ勝局のみを抽出したデータを、次に敗局も含めたデータをrstanを用いて分析した。 研究仮説は「第2群の母平均の方が第1群の母平均より大きい」としておこう。

勝局のみの分析

上記ソースから抽出した一致率を羅列すると以下のようになる。

  • 6月以前
    78.9%, 53.8%, 52.8%, 54.5%, 67.6%, 68.2%, 63.0%

  • 7月以降
    74.1%, 80.6%, 83.3%, 80.8%, 81.8%

まずは一致率が正規分布に基づくと仮定してモデルを立てる。

library(rstan)
miura <- "
data {
  int<lower=0> n1;                        #6月以前のデータ数
  int<lower=0> n2;                        #7月以降のデータ数
  real<lower=0> x1[n1];                  #6月以前のデータ
  real<lower=0> x2[n2];                  #7月以降のデータ
  real mL; real mH; real sL; real sH;  #事前分布
}
parameters {
  real<lower=mL,upper=mH> mu1;        #平均(範囲指定)
  real<lower=mL,upper=mH> mu2;        
  real<lower=sL,upper=sH> sigma1;       #標準偏差(範囲指定)
  real<lower=sL,upper=sH> sigma2;       
}
model {
  x1 ~ normal(mu1,sigma1);                #正規分布に従うと仮定
  x2 ~ normal(mu2,sigma2);                
}
"

   
データとパラメータを設定してMCMCを実行。

n1 <- 7
n2 <- 5
x1 <- c(78.9,53.8,52.8,54.5,67.6,68.2,63.0)
x2 <- c(74.1,80.6,83.3,80.8,81.8)
mL <- 0
mH <- 100 #平均の範囲は0から100とする
sL <- 0
sH <- 100 #標準偏差の範囲は0から100とする
prob <- c(0.025, 0.05, 0.5, 0.95, 0.975)
see <- 1234
cha <- 3
war <- 1000
ite <- 21000
s3 <- stan_model(model_code=miura)
dat <- list(n1=n1, n2=n2, x1=x1, x2=x2)
fit <- sampling(s3, data=dat, seed=see, chains=cha, warmup=war, iter=ite)

 
得られたリストから平均のデータを抽出、研究仮説「第2群の母平均の方が第1群の母平均より大きい」が正しい確率を計算してみる。

mu1list <- extract(fit, "mu1")
mu2list <- extract(fit, "mu2")
difference <- unlist(mu2list[[1]])-unlist(mu1list[[1]])
whether <- ifelse(difference > 0, 1, 0)
mean(whether)

 
結果は以下のようになる。

> mean(whether)
[1] 0.9947667

 
よって、勝局のみを考えた場合、研究仮説「第2群の母平均の方が第1群の母平均より大きい」が正しい確率は99.5%となる。

敗局も含めた分析

敗局も含めた一致率の羅列は以下のようになる。

  • 6月以前
    78.9%, 53.8%, 52.8%, 54.5%, 67.6%, 68.2%, 63.0%, 40.0%, 36.0%

  • 7月以降
    74.1%, 80.6%, 83.3%, 80.8%, 81.8%, 58.3%, 48.5%, 48.6%

n1. <- 9
n2. <- 8
x1. <- c(78.9,53.8,52.8,54.5,67.6,68.2,63.0,40.0,36.0)
x2. <- c(74.1,80.6,83.3,80.8,81.8,58.3,48.5,48.6)

dat. <- list(n1=n1., n2=n2., x1=x1., x2=x2.)
fit. <- sampling(s3, data=dat., seed=see, chains=cha, warmup=war, iter=ite)
mu1list. <- extract(fit., "mu1")
mu2list. <- extract(fit., "mu2")

difference. <- unlist(mu2list.[[1]])-unlist(mu1list.[[1]])
whether. <- ifelse(difference. > 0, 1, 0)
mean(whether.)

 
結果は以下のようになる。

> mean(whether.)
[1] 0.9271167

 
よって、敗局も考慮した場合、研究仮説「第2群の母平均の方が第1群の母平均より大きい」が正しい確率は92.7%と求められる。

まとめ

以上より、6月以前の一致率の平均値が以降の平均値より大きい確率は

  • 勝局のみを対象とした場合:99.5%

  • 敗局も含んだ場合:92.7%

となった。

私はこの数字に関して一定の閾値を上回る或いは下回ると議論することはしないし、一将棋ファンである私がこれに基づいて何ら見解を表明する立場にないことは明らかである。無責任にデータを放り投げて「解釈は自由」とだけ述べ、ひとまず筆をおくこととしよう。