TopCoder SRM420 Div2 250Pts

このTopCoderの問題はこちらで見ることができる(要TopCoder登録 & 問題文は英語)。問題文についておおまかに説明する。

カードのデッキがある。それぞれのカードは黒ないしは赤である。以後、Bという文字を黒のカード、Rという文字を赤のカードとする。デッキの上から下に向かってカードを並べることによって、deckという文字列を生成する(訳注：つまり、deckの中身はRとBという文字列を利用して構成する)。次に、deckのすべてのカードを挿入して新しいデッキを作成したい。まず始めに、新しいデッキは空であるとする。次に古いデッキの上から下に向かってカードを処理していく。deckのi文字目に対応するカードを処理するというのは、古いデッキからカードを取り出し、新しいデッキにabove[i]だけ上にカードがある位置にカードを挿入するということになる。

deckという古いデッキと、above[]という新しいデッキへの挿入位置が与えられたとき、deckと同様な形式の文字列に変換された新しいデッキを表す文字列を返すメソッドを作成せよ。

例えばdeck="BRBRR"、above={0, 0, 1, 0, 3}であれば、"RRBRB"となる。これは初めにBを挿入(状態：B)、次に上に0枚がある位置にRを挿入(状態：RB)、次に上に1枚ある位置にBを挿入(状態：RBB)、...とすることによって得られる結果である。

私の実装はこちら。

public class DeckRearranging {

 public String rearrange(String deck, int[] above) {
  String ret = "";
  for( int i=0 ; i<above.length ; i++ ){
   String first = ret.substring(0, above[i]); // 新しいdeckを分けたときに上にあるカード群
   String second = ret.substring(above[i]); // 新しいdeckを分けたときに下にあるカード群
   ret = first + deck.charAt(i) + second; // 新しいdeckにカードを挟み込む
  }
  return ret;
 }

}

得点は221.63/250、1回のsubmitでシステムテストクリア。中央値は約212点。StringBufferのinsertを利用するのもいいかもしれません。

ニュートン法による極値計算

テキスト「これなら分かる最適化数学」の86ページにあるニュートン法をPythonで実装した。通常のニュートン法は関数f(x)に対し、f(x)=0なるxを求める反復計算を行うものであるが、反復の式を少し変更するだけで、極値(f'(x)=0なるx)を求めることができる。詳しい理論はテキストに譲るが、x <- x - f'(x)/f''(x)とし、収束するまで繰り返せば、極値にたどり着くことができる。

ニュートン法で極値を求めるにあたり、注意しなければならないのは、コードの更新式からも見当がつくが、途中で2回微分が0にならないようにすることである。

以下、f(x)=x^3-2*x^2+x+3に対して、極値が求まるか試してみた。なお、f(x)が極値をとるxの値は解析的に解くこともでき、x=1、1/3の2つとなる。

# -*- encoding: utf-8 -*-

'''
ニュートン法による極値の計算
2回微分が存在する場合、勾配法よりも効率よく極値を計算できる
'''

from sympy import * # 微分操作のために必要
import itertools as it # mapで複数の引数を与えたいので利用する

def newton(init, vari, func, delta=10**-6):
    '''
    funcは極値を探す対象の関数、variは微分する変数、deltaは収束判定のための定数
    返す値は見つかった極値の1つである
    '''
    df = diff(func, vari)
    d2f = diff(df, vari)
    x = init
    while True:
        xvar = x
        # 2回微分の微分係数が0にならないように注意が必要
        x = xvar - float(limit(df, vari, xvar)) / float(limit(d2f, vari, xvar))
        if abs(x - xvar) < delta:
            break
    return round(x, 5) # deltaの値によって有効桁数が変わるが、簡単のため5で決めうち

x = Symbol('x')
f = x**3 - 2*x**2 + x + 3
start_num = range(-5, 5) # 初期値を複数設定しておく
print map(newton, start_num, it.repeat(x, len(start_num)),
          it.repeat(f, len(start_num)))

最後のmapを利用することで、一度に複数の結果を実行、表示することができるようになる。mapに与える引数については、繰り返し実行できる形式にしておく必要がある。

実行結果は、[0.33333, 0.33333, 0.33333, 0.33333, 0.33333, 0.33333, 1.0, 1.0, 1.0, 1.0]となる。それぞれ初期値を-5、-4、...、4に変更しながら実行した結果である。すべて1/3と1のうち、近いほうの極値に収束していることが確認できる。

勾配法を用いた最適化計算

テキスト「これなら分かる最適化数学」の81ページにある勾配法(gradient method)アルゴリズムをpythonで実装した。勾配法とは、関数f(x)を最大にするxの値に近いと思われる値を初期値とし、必ず関数値f(x)が増加するようにし、かつステップ幅をなるべく大きくなるようその幅を調整することにより、できるだけ高速に関数を最大にするxを求めるというものである。

なお、今回扱う勾配法においては、変数の値がとりうる範囲について、関数の微分係数が0になるのは1箇所であると仮定する(複数極値があると、どこに向かうか分からないため)。また、微分できない箇所がある関数では不都合が生じるので、これも考えないことにする(例えばf(x)=|x|のx=0の点のように、尖っている状況下などがある)。

関数の最大値を求めるにあたっては、微分係数が0になるという性質を利用している。このため、微分という操作を行う必要であるが、この問題はsympyモジュールを利用することによって解決できる。sympyモジュールはこちらからダウンロードできる。なお、記事執筆時点のsympyのバージョンは0.7.1である。

以下、「2次関数f(x)」と「正規分布の確率変数と関係がある部分のみ取り出した関数g(x)」の2つの例を用いて、勾配法で正しく関数を最大化するようなxを推定できるか試してみた。結果としては、(制約が強いこともあるが)ほぼ正確に推定できているということになった。

# -*- encoding: utf-8 -*-

from sympy import *

# 符号関数
def sign(value):
    if value > 0: return 1
    if value < 0: return -1
    return 0

# 勾配法を計算するメインルーチン
# initは初期推定値、funcは最大値を求めたい関数、stepは初期移動量、
# epsは収束した(ほぼf'(x)が0となった)と判定する閾値である。
# 返り値は最大値となるxの値に相当する。
def gradient_method(init, func, step = 10**-7, eps = 10**-6):
    # diffは微分の関数、これで第一引数の関数をxで微分したことになる
    # 事前にxをSymbol関数でシンボル化しておく必要がある
    df = diff(func, x)
    xx, h = init, step
    while True:
        h = sign(limit(df, x, xx)) * abs(h)
        X, Xdash = xx, xx + h
        # 関数に値を入れるという方法が見つからないので
        # 極限の値を求めている
        if limit(func, x, X) < limit(func, x, Xdash):
            while limit(func, x, X) < limit(func, x, Xdash):
                h *= 2.0
                X, Xdash = Xdash, X + h
            xx, h = X, h / 2.0
        else:
            while limit(func, x, X) > limit(func, x, Xdash):
                h /= 2.0
                Xdash -= h
            xx, h = Xdash, 2.0 * h
        if abs(limit(df, x, xx)) <= eps:
            break
    return xx

x = Symbol('x') # xを記号として扱う
f = -x**2 + 2*x + 3
g = exp(-(x-2)**2)
print gradient_method(0, f) # 1.0000002、なお正解は1である
print gradient_method(0.5, g) # 2.0000002、なお正解は2である

TopCoder SRM419 Div2 250Pts

このTopCoderの問題はこちらで見ることができる(要TopCoder登録 & 問題文は英語)。問題文についておおまかに説明する。

柱状図(訳注：棒グラフと思って差し支え無いです)とは水平方向にn個の柱が配置されたものである。各柱は幅が1であり、各柱はx軸にくっついている。柱の高さはa[]という整数型の配列で与えられる。a[i]はi番目の柱の高さを表している(訳注：より詳しくはリンク先の画像をご確認ください。隣合う柱はくっついている図です)。与えられた柱状図の周囲の長さを返すメソッドを作成せよ。なお、a[i]は1以上50以下の整数になるものとする。例えば、aが{3, 2, 1}であれば、x軸に水平な方向の周囲の長さは柱の上側、下側を合わせて3*2=6、側面の長さは左から順に3+1+1+1で合計6なので、それらを合わせて12になる。

私の解答はこちら。

public class ColumnDiagramPerimeter {

 public int getPerimiter(int[] a) {
  int perimiter = 0;
  int prev = 0;
  for( int i=0 ; i<a.length ; i++ ){
   perimiter += Math.abs(a[i] - prev);
   prev = a[i];
  }
  perimiter += (a[a.length-1] + a.length * 2);
  return perimiter;
 }

}

得点は207.19、2回目のsubmitでシステムテストクリア。途中で提出のページに接続できなくなり、提出が遅れたのが残念。

TopCoder SRM418 Div2 250Pts

このTopCoderの問題はこちらで見ることができる(要TopCoder登録 & 問題文は英語)。問題文についておおまかに説明する。

あなたは真剣に戦略ゲームをするプレーヤーなので、敵の防御塔を攻撃するという、もっとも一般的な問題を解くことにした。攻撃をする前にはmyUnitsという数の兵士を得る。各兵士は1回のラウンドにおいて、1つの塔に1のダメージを与えることができる。敵は兵士を持たない。しかしながら、体力がhpTの塔をnumT所持している。各塔は1回のラウンドにおいて、attackTだけの兵士を殺すことができる。1つのラウンドは次の順で行われる。

1. 自軍の各兵士は塔を1つ選び、塔に1のダメージを与える。塔がすべての体力を失うと、塔は破壊される。各兵士ごとに独立して攻撃する塔を選べるものとする。
2. 敵の攻撃になる。敵はk*attackTだけの自軍の兵士を殺すことになる。ここで、kは残っている塔の数である。

あなたの課題はすべての塔を破壊することである。もしそれが可能であれば、それを実行するのに必要な最小限のラウンド数を返せ。もし不可能ならば-1を返すようなメソッドを作成せよ。

私の解答はこちら。

public class Towers {

 public int attack(int myUnits, int hpT, int attackT, int numT) {
  int nTurn = 1; // 攻防を行った回数
  int totalHP = numT * hpT; // 塔の合計体力
  int remainT = numT; // 残っている塔の数
  int remainU = myUnits; // 残っている自軍の兵士数
  while( true ){
   totalHP -= remainU; // 塔の体力を減らす
   if( totalHP <= 0 ) return nTurn; // 塔の体力の合計が0以下なら全部破壊されたことになる
   remainT = (int)Math.ceil(totalHP*1.0/hpT); // 残っている塔の数を計算
   remainU -= remainT * attackT; // 自軍の兵士の数を減らす
   if( remainU <= 0 ) return -1; // もし自軍の兵士が全滅したら-1を返す
   nTurn++;
  }
 }

}

得点は196.74/250、1回のsubmitでシステムテストクリア。中央値は約135点。すべての塔を破壊するにはできるだけ塔を集中的に攻撃するというのが最善になることを利用します。タワーの残り体力を計算する際に、totalHPに1.0を掛けていなかったことに気付くのに少し時間がかかりました。

TopCoder SRM417 Div2 250Pts

このTopCoderの問題はこちらで見ることができる(要TopCoder登録 & 問題文は英語)。問題文について意訳で説明する。

与えられたXという整数に対し、Rev(X)というXの桁を反転させた整数を得ることができる。例えばRev(123)=321、Rev(100)=1である。xとyという二つの正の整数が与えられたとき、Rev(Rev(x) + Rev(y))を返すメソッドを作成せよ。

私の解答はこちら。

public class ReversedSum {

 public int getReversedSum(int x, int y) {
  StringBuffer xx = new StringBuffer("" + x);
  StringBuffer yy = new StringBuffer("" + y);
  int s = Integer.parseInt(xx.reverse().toString()) +
                        Integer.parseInt(yy.reverse().toString());
  StringBuffer ss = new StringBuffer("" + s);
  return Integer.parseInt(ss.reverse().toString());
 }

}

得点は243.24/250、1回のsubmitでシステムテストクリア。中央値は約233点。エレガントさを求めるのであればRev(X)という関数を定義するのが良さそうです。

TopCoder SRM416 Div2 250Pts

このTopCoderの問題はこちらで見ることができる(要TopCoder登録 & 問題文は英語)。問題文について意訳で説明する。

英語の典型的なテキストについて、いくつかの文字が他の文字よりも頻繁に現れるということが一般的に知られている。例えば、英語の長文においては、約12.31%の文字がeになる。言語学部の友人が、与えられたテキストの中で最も頻繁に現れる文字が知りたいということで、手伝ってくれるよう頼んできた。テキストはtext[]という文字列型の配列で与えられる。このとき、テキストでもっとも頻繁に現れる文字を返すメソッドを作成せよ。もし一番頻繁に現れる文字が複数ある場合、それらをアルファベット順で返すようにせよ。なお、文字列は小文字のアルファベットと空白しかないものとする。

私の解答はこちら。

public class MostCommonLetters {

 public String listMostCommon(String[] text) {
  int[] times = new int[26];
  for( int i=0 ; i<text.length ; i++ ){ // ここで文字の頻度集計
   for( int j=0 ; j<text[i].length() ; j++ ){
    char c = text[i].charAt(j);
    if( c != ' ' ){
     times[c-'a']++;
    }
   }
  }
  int maxTimes = -1;
  for( int i=0 ; i<times.length ; i++ ){ // 最大出現回数を調べる
   maxTimes = Math.max(maxTimes, times[i]);
  }
  String ret = "";
  for( int i=0 ; i<times.length ; i++ ){ // 出現回数が最大となる文字を集める
   if( times[i] == maxTimes ){
    char c = (char)('a'+i);
    ret += c;
   }
  }
  return ret;
 }

}

得点は241.25/250、1回のsubmitでシステムテストクリア。中央値は約186点。素直な集計の問題です。

TopCoder SRM415 Div2 250Pts

このTopCoderの問題はこちらで見ることができる(要TopCoder登録 & 問題文は英語)。問題文について意訳で説明する。

あなたの趣味は郵便消印を集めることである。N個の異なる消印があり、0からN-1まで番号が付けられている。それらの価値はprices[]という整数型の配列で与えられており、添字iが消印iの価値を表すものとする。あなたの目標はできるだけ異なる種類の消印を集めることである。既に持っている消印はhave[]という整数型配列で与えられている。また、初期状態では所持金0の状態とする。消印を売り、お金を得て異なる消印を買うということができる。集められる消印の種類数の最大値を返すメソッドを作成せよ。

私の解答はこちら。

import java.util.Arrays;

public class CollectingUsualPostmarks {

 public int numberOfPostmarks(int[] prices, int[] have) {
  int value = 0;
  for( int i=0 ; i<have.length ; i++ ){
   value += prices[have[i]];
  }
  Arrays.sort(prices);
  int count = 0;
  int current = 0;
  for( int i=0 ; i<prices.length ; i++ ){
   current += prices[i];
   if( current > value ) break;
   count++;
  }
  return count;
 }

}

得点は191.82/250、問題文を読み替えて、全部売って安い方から買っていくと問題を置き換える発想の転換すると非常に見通しが良くなります。

一様分布から標準正規乱数を発生させる

標準正規分布に従う乱数を発生させる簡単な方法として、一様乱数をいくつか発生させてその和をとり、期待値で引き算するという方法がある。そこで、得られた値が本当に標準正規分布といって差し支えないか調べてみることにした。以下のコードでは、6つの一様乱数(取りうる値の範囲は[0,1])とし、それらの和の期待値3を引いたものを大量に発生させ、値の分布が一様乱数であるか検定したものである。6としているのは、とりうる値を[-3, 3]とするためである。標準正規分布であれば、多くの場合、そこから発生する乱数は[-3, 3]の範囲に収まるためである。

# 一様乱数から正規乱数を生成するプログラム

set.seed(1)
pts <- 6 # 一つの標準正規乱数を発生させるのに用いる一様乱数の数
N <- 5000 # 発生させる正規乱数
num.data <- N*pts
rn <- runif(num.data)
rnd.mat <- matrix(rn,nrow=N) # 一様乱数の行列
num <- rowSums(rnd.mat)-pts/2

hist(num,breaks=20,xlim=c(-4,4),xlab="",ylab="freq",main="")
oldpar <- par(new=T)
curve(dnorm,-4,4,xlab="",ylab="",axes=FALSE)
par <- oldpar

# コルモゴロフ・スミルノフ(KS)検定
# 標準正規分布 N(0, 1)とデータを比較している
ks.test(num, "pnorm", mean=0, sd=1)

コルモゴロフ・スミルノフ検定は、データがある分布に従っているかを検定するというものである。上のコードは1標本KS検定に相当する。帰無仮説はデータは分布から発生したものではない、対立仮説はデータは分布から発生したものとなる。統計量は従うと想定される分布の累積分布関数と、データのそれを比較した結果から得られる差の値に基づく。統計量が小さくなると、帰無仮説が棄却されるということになる。結論は、KS検定ではp値はほぼ0となり、帰無仮説は標準正規分布ではないということになった。ただし、中心極限定理からデータ数を6から増やせば、標準正規分布に従うことが想定される。そこで、pts<-12としてみると標準正規乱数といえるようなp値となった（有意水準5%で帰無仮説を棄却できる）。必要なptsのサイズは標準化する分布(今回は一様分布)の対称の程度に依存していることに注意。今回利用した一様分布U(0,1)は0.5で対称なので、比較的少ないデータ数で標準正規分布が生成できるのである。

TopCoder SRM414 Div2 250Pts

このTopCoderの問題はこちらで見ることができる(要TopCoder登録 & 問題文は英語)。問題文について意訳で説明する。

あなたはレストランのサービス部門のマネージャである。夕方のコースでは、顧客のグループがやってきて、その人たちをテーブルに配置しなければならない。レストランのテーブルは異なった人数の客に対応するために、さまざまなサイズのテーブルを用意している。テーブルへの客の配置は次のような手順で行う。客が現れたとき、まだ客で埋まっていないテーブルの内、そのグループ以上の人数を収められる最小のテーブルを客のために確保する。もし適切なテーブルが複数あれば、その中から任意に一つ割り当てる。もし適切なテーブルが一つもない場合は、顧客は即座に踵を返してその日の夕方は二度と現れない。グループが食事を終え、離れるとき、テーブルは利用できるようになり、新しい顧客を配置することができるようになる。

あなたはこの方法でテーブルに客を割り当てたとき、何人の顧客に逃げられるか知りたい。テーブルのサイズはtable[]という整数型配列で、各要素は収容可能な人数を表す。グループの人数、グループの到着時刻、出発時刻はそれぞれgroupSizes[i]、arrivals[i]、departures[i]である。i番目の要素はすべてグループiに関する情報である。グループは到着時刻順に並べられており、どのグループも同時に来ることはないものとする。もし、あるグループの出発時刻と別のグループの到着時刻が一致した場合は、出発するグループのテーブルは利用可能になっているものとする。逃げられた顧客の総数を返すメソッドを作成せよ。

私の解答はこちら。

import java.util.Arrays;

public class RestaurantManager {

 public int allocateTables(int[] tables, int[] groupSizes, int[] arrivals, int[] departures) {
  Arrays.sort(tables);
  boolean[] isOccupied = new boolean[tables.length];
  // テーブルにいるグループの出発時刻。
  // i番目の要素がテーブルiを表す。値0はそのテーブルに誰も座ってないことを意味する
  int[] depTime = new int[tables.length];
  int nTurnedAway = 0;
  for( int i=0 ; i<arrivals.length ; i++ ){
   for( int j=0 ; j<tables.length ; j++ ){ // 来た人の時刻により、埋まっていたテーブルを解放
    if( arrivals[i] >= depTime[j] && depTime[j] > 0 ){
     isOccupied[j] = false;
     depTime[j] = 0;
    }
   }
   boolean caught = false;
   for( int j=0 ; j<tables.length ; j++ ){ // 空いているテーブルの探索
    // 最初に見つかったところに顧客を割り当てる(事前にテーブルサイズでソートしておくのがキモ)
    if( isOccupied[j] == false && tables[j] >= groupSizes[i] ){
     isOccupied[j] = true;
     depTime[j] = departures[i];
     caught = true;
     break;
    }
   }
   if( caught == false ){
    nTurnedAway += groupSizes[i];
   }
  }
  return nTurnedAway;
 }

}

得点は201.97/250、1回のsubmitでシステムテストクリア。中央値は約125点。テーブルを小さい順に並び替えておくことで、空いている必要最上限のテーブルサイズを探索するのが容易になります。

TopCoder SRM413 Div2 250Pts

このTopCoderの問題はこちらで見ることができる(要TopCoder登録 & 問題文は英語)。問題文について意訳で説明する。

地下鉄の電車は人々をある駅から隣の駅へ高速に動かすことができる。二つの隣り合う駅について、lengthメートルであることが知られている。安全のため、電車はmaxVelocityメートル/秒よりも速く移動することは許されない。また、快適さのため、加速度の絶対値はmaxAccelerationメートル/秒^2を超えてはならないものとする。電車は0メートル/秒から動きだし、隣の駅で停車する(つまり駅に到着したときに0メートル/秒の速度になる)。隣の駅へ到着するのに必要な時間の最小値を返すようなメソッドを作成せよ。

私の解答はこちら。

public class Subway2 {

 public double minTime(int length, int maxAcceleration, int maxVelocity) {
  double toMax = maxVelocity * 1.0 / maxAcceleration;
  double s = maxVelocity * toMax;
  if( length < s ){
   double half = length / 2.0;
   return Math.sqrt(half * 2 / maxAcceleration)*2;
  }else{
   return toMax * 2 + (length - s) / maxVelocity;
  }
 }

}

209.07/250、1回のsubmitでシステムテストクリア。高校の物理学に出てくる運動の問題ですね。最高速に達する(コード中の分岐では後半)か、達しない(コード中の分岐では前半)かで場合分けすると簡単です。

2012年7月の学習記録

2012年7月に学習した主な本一覧。

入門・演習 数理統計(pp.68～97、65問)

積率母関数や条件付き期待値の問題などを解く。

マスタリングTCP/IP(pp.231～338読了)

この本の範囲はエンジニアなら知っておいて損はないように思いました。とりあえず自分の周囲半径3メートル程度で発生したネットワークトラブルは自分の力で原因究明できるようになりたいものです。

テキストデータの統計科学入門(pp.13～134)

テキストをうまく指標に落として検定等に持ち込む方法について学習中。個人的にはRでネットワークを作るパッケージの紹介が参考になりました。

統計解析入門（pp.130～148)

信頼区間は区間が最小になるように選んでいると知ってナルホドと思いました。

入門自然言語処理(pp.62～110)

もう少し深く読みたいです。個人的には演習問題の解答を一通り作ってしまいたい。

初めてのPython(第13～16、22～23章)

どちらかというとリファレンス的な本なので、実践的なテキストを読んだ方が楽しそう。

その他、実践ビジネス英語、入門ソーシャルデータなどを少々掻い摘む。