それが僕には楽しかったんです。

僕と MySQL と時々 MariaDB

Newton-Raphson法と単位円内接正多角形を使った円周率近似をScalaでやってみた

やってみた

C言語縛りは面白くないので例によってScalaでやってみた

こんなやつをコードにした

まずNewton-Raphson法は次の式におけるxの値を求める

 \displaystyle
f(x) = g(x) - C = 0

このあと登場するけど、単位円内接正多角形を使う近似では漸化式をつかっていくわけだけども
その中で平方根を計算しないといけないのでCを任意実数定数として関数g(x)を次のようにする

 \displaystyle
g(x) = x^2

すると次の式におけるxを近似的に求めることになる

 \displaystyle
f(x) = x^2 - C = 0 \\
x^2 - C = 0 \\
x^2 = C

Newton-Raphson法ではこれを次の式を満たす様に微分を使いながら近似していく

{ \displaystyle
x_{n+1} = x_n - \frac{f(x_n)}{\frac{df}{dx}(x_n)}
}

ここまでNewton-Raphson法

次に単位円を用いて円周率を求める。
しかし、円周がわからない。
これは円周率が不明であるためである。

そこで単位円内に内接正多角形を作って円周とみなす。
これは漸化式で表現でき正n角形のnが大きくなればなるほどその周の長さが円周の長さに近づく。

その内接正多角形の一辺の長さが次の式となる
{ \displaystyle
L_{n+1} = \sqrt{ 2-2\sqrt{1-(\frac{L_n}{2})^2} } \: or \:  L_{n+1} = \sqrt{2-\sqrt{4-(L_n)^2}}
}

というわけでコードにしてみた

実際に書いたコード

object calcPI{

    val max = 20

    def main(args:Array[String]):Unit = {
        var Ln = newtonRaphson(2)//L_1

        var n = 0
        var pow = 2
        while( n < max ){
            println( pow * Ln )
            Ln = newtonRaphson(2 - newtonRaphson(4-Ln*Ln))
            n += 1
            pow *= 2
        }
    }

    def newtonRaphson(c:BigDecimal):BigDecimal = {
        var x:BigDecimal = c
        var dx:BigDecimal = -func(x,c) / dfunc(x)

        while( 1.0e-16 < dx.abs  ){
            x += dx
            dx = -func(x,c) / dfunc(x)
        }

        return x
    }

    def func(x:BigDecimal, c:BigDecimal):BigDecimal = x*x - c
    def dfunc(x:BigDecimal):BigDecimal = 2.0*x
}

これを実行して、出力結果をdatファイルにおとしてgnuplotでグラフにするとこんな感じになる

f:id:RabbitFoot141:20170504120319p:plain

BigDecimal使ってはいるけど精度的にはまだまだかなと思うけどある程度できた。

Newton-Raphson法を用いた円周率近似

Newton-Paphson法の基礎原理

Newton-Raphson法とは次の式を満たすxの値を見つけるアルゴリズム
{ \displaystyle
f(x) = g(x) - C = 0
}
ここではg(x)を任意の微分可能な関数としてCを実数定数とする。
つまりNewton-Rapshon法は次の式を満たすxの値を見つけるアルゴリズムとなる。

{ \displaystyle
g(x) = C
}

つまり以下の様にg,Cが与えられた時は2の平方根を求めることになる。

{ \displaystyle
g(x) = x^2, C = 2
}

実際の計算

基礎原理はうえで紹介したので実際に解を近似して求める。というのも上記の式を満たすxがすぐに正確に見つかるわけではないことによる。
そこで与えられたf(x)を一回微分して接線の方程式を求め、その接線とx-軸の接点におけるx座標で再びf(x)に戻り微分し接点のx座標を求める…。
これを再帰的に繰り返すことで求めたいxに近似していくこととなる。
その漸化式は次のようになる。

{ \displaystyle
x_{n+1} = x_n - \frac{f(x_n)}{\frac{df}{dx}(x_n)}
}

この漸化式をC言語を用いたコードとして表現していく。
ここでは平方根を求めることにする。

コード化する

ひとまず次のプログラムをベースにする。

#include<stdio.h>

int main(){
    return 0;
}

ここで関数fとそれを一回微分したdfを関数として作る。

#include<stdio.h>
#include<stdlib.h>
#include<math.h>

//プロトタイプ宣言
double Func(double x,double C);//変数xとCが必要となる
double dFunc(double x);//fを微分しているのでxのみ必要

int main(){
    return 0;
}

double Func(double x, double C){
    return x*x - C; //x^2-Cを計算して返す
}

double dFunc(double x){
    return 2.0*x; //2xを返す
}

これで準備は出来た。

次に微分して…。となるがCで微分を扱うのは難しいので数学的に考えてfとxのみの式に直す。
それは以下になる。

{ \displaystyle
\frac{df}{dx}(x_n) = \frac{f(x_n)}{x_n-x_{n+1}}
}

というわけでNewton-Raphson法を用いて平方根を求める。ここでは誤差が{10e-6}になるまで続けることにする。

#include<stdio.h>
#include<stdlib.h>
#include<math.h>

//追加した
#define EPS 10e-6

//プロトタイプ宣言
double Func(double x,double C);//変数xとCが必要となる
double dFunc(double x);//fを微分しているのでxのみ必要

double Newton(double C);//平方根を求めたい値を実数で受け取る

int main(){

    printf("%12.10f\n",Newton(2) );

    return 0;
}

//Newton-Rapshon法を用いてcの平方根を求める。
double Newton(double C){
    
    double x = C;
    double dx = -Func(x,C)/dFunc(x);
    
    while( fabs(dx) > EPS){
        x += dx;
        printf("%12.10f\n",x);
        dx = -Func(x,C)/dFunc(x);
    }
    return x;
}

double Func(double x, double C){
    return x*x - C; //x^2-Cを計算して返す
}

double dFunc(double x){
    return 2.0*x; //2xを返す
}


これで平方根をNewton-Raphson法を使って求めることができる。
ただし、誤差があることに注意。

円周率近似

円周率近似も漸化式を用いて行う。
http://oto-suu.seesaa.net/article/291513400.html
上記のサイトがそれを紹介しているが半径が{\frac{1}{2}}なので、半径を単位円において考えることにする。
基本的には正n角形の周を求め円周に近似していくことになる。

その場合、漸化式は次のようになる。
{ \displaystyle
L_{n+1} = \sqrt{ 2-2\sqrt{1-(\frac{L_n}{2})^2} } or L_{n+1} = \sqrt{2-\sqrt{4-(L_n)^2}}
}

ソースコードにする

これをソースコードにするとこうなる

#include<stdio.h>
#include<stdlib.h>
#include<math.h>

#define EPS 1.0e-8
#define POW 2

double Newton(double);
double Func(double, double);
double dFunc(double);

int main(){

    double Ln = Newton(2.0);//L_1 = sqrt(2)
    printf("%12.10f\n",Ln);
    int n = 1;
    int pow_base = 2;
    int n_max = 20;

    while( n < n_max ){ 
        
        printf( "%12.10f\n", pow_base * Ln );
        Ln = Newton( 2-Newton(4 - pow(Ln, 2)) );
        n += 1;
        pow_base *= POW;

    }

    return 0;
}


//平方根を求める
double Newton(double C){
    
    double x = C;
    double dx = -Func(x, C) / dFunc(x);
    
    while( fabs(dx) > EPS ){
        x += dx;
        dx = -Func(x, C) / dFunc(x);
    }

    return x;
}


double Func(double x, double C){
    return x*x-C;
}

double dFunc(double x){
    return 2.0 * x;
}

これで円周率の近似値が求められた。

githubでコミットしても芝生が生えなかったわけ

昨日、githubの自分のレポジトリにコミットしてたら唐突に芝生が生えなくなった。

それで色々ググった結果

  • masterブランチじゃないとだめ
  • ローカルレポジトリの登録アドレスがgithubアカウントのメールアドレスになっていないとだめ
  • forkしたやつもだめ

みたいなあれこれが書いてあって、それ以外の原因は見つからなかった。

がしかし、自分の環境ではmasterブランチにコミットしていたし、メアドも合ってるし、自分で作成したレポジトリだしってことで
沼に陥ったので、github supportにこれらの問題をぶん投げた。

すると返信があって

「ラグってるだけだと思うから、24時間以上して芝生生えなかったらまたメールしてね」

と来たので今朝見たら、本当にラグってただけだった。

なんか芝生もすぐに反映されなくて、数時間ラグることもあるらしい。

遺伝的アルゴリズムとその応用

遺伝的アルゴリズム ~Genetic Algorithm~

遺伝的アルゴリズム(Genetic Algorithm)とは、近似解を探索するアルゴリズム
主にランダム性を持って要素を変化させる。進化的アルゴリズムの一つでもある。

このアルゴリズムはデータを遺伝子として表現し、一つ一つを個体とする。
その個体を複数用意して、適応度計算をし低いものを交叉や組み換え、突然変異させ解を探っていく。
個体を何世代も変化させていくので、組み合わせ最適化問題(後述)などで効果を発揮する。

遺伝的アルゴリズムの基本原理

遺伝的アルゴリズムの基本原理を説明する。

ここでは n \in \mathbb{N} 個の個体と g \in \mathbb{N}世代で行うと仮定する。

  • 初めに現世代と次世代の個体が入る集合をそれぞれ一つずつ用意しておく。
  • 現世代にn個の個体をランダムに生成する。
  • 評価関数を使用し、現世代の個体が持つ適応度を計算する
  • ある確率において、次の3つのどれかを行い結果を次世代に保存する。これを次世代個体がn個になるまで行う。
    • 個体を2つ選択し交叉させる
    • 個体を一つ選択し、突然変異させる
    • 個体を一つ選択し、そのままコピーする
  • 次世代個体を現世代にコピーし評価し直す
  • これまでの手順をg回繰り返した後、最も適応度の高い個体が解になる確率が高い。

遺伝的アルゴリズムにおける個体操作

これまで出てきたのは次の3つ

  • 選択
  • 交叉
  • 突然変異

選択

選択は自然淘汰をモデル化していて、適応度に基づいて個体を削除したり作成する操作のこと。
この選択の手法として次の3つが挙げられる

  • ルーレット選択
  • ランキング選択
  • トーナメント選択

これらが代表的な選択だが、ここに出てきた選択を使用せず
高い適応度の個体を一定数そのまま残すという手法を取ることでも似た結果が得られる。
ただし、選択アルゴリズムを使うと初期収束の原因になったり
高い適応度の個体を残し過ぎると解の多様性が失われるなどの問題がある。

交叉

交叉、つまり組み換えは個体のもつ情報を組み替え適応度に変化をもたせるもの。
選択では適応度自体が変化する訳ではないので
遺伝的アルゴリズムにおいて交叉は最適解に近づけるための重要な遺伝操作といえる。
この交叉というのにも以下のような種類がある。

  • 一点交叉
  • 二点交叉
  • 多点交叉
  • 一様交叉

交叉は重要な操作であるから、この中から最適なものを選ぶ必要がある。
一番実装しやすいのは1/2の確率で要素を入れ替える一様交叉だと思う。

突然変異

突然変異は選択と同様に、生物に見られる遺伝的な特徴をモデルにしたもの。
選択は生物の自然淘汰がモデルだったのに比べ、突然変異はその名の通り突然変異をモデルにしたものである。

突然変異は交叉や選択に関わらず一定の確率で行う必要があるがその確率は0.1%~1%台が望ましく
高くても数%以内にするべきと言われている。
この確率が高すぎると、遺伝的アルゴリズムではなくちょっと面倒な乱択アルゴリズムになってしまい、最適解に収束しにくくなる。

遺伝的アルゴリズムの問題点

遺伝的アルゴリズムには代表的な問題点が存在する。

  • 初期収束
  • ヒッチハイキング

である。

初期収束

これは最初の方の世代で偶然適応度がものすごく高い個体が生成され。
世代が進むにつれて、その個体が現世代次世代を通して爆発的に増加してしまうこと。
中途半端な段階で収束してしまうと、最適解を得ることが難しくなる。
この時は突然変異の確率や、選択の確率を見直す必要がある。

ヒッチハイキング

交叉において、最適解の生成を妨げる遺伝をしてしまう現象。
これは等確率な交叉を行う一様交叉を使用することで防ぐことができる。

応用例

遺伝的アルゴリズムを使用して、「ナップサック問題」を解いてみる。
これは組み合わせ最適化問題の一つである。

問題

上限が6kgのナップサックに対して値段と重量が以下の品物を入れる。
そのとき、上限重量以下で値段が最高になる組み合わせを求める。

A:1kg 100-yen
B:2kg 300-yen
C:3kg 350-yen
D:4kg 500-yen
E:5kg 650-yen

余談だが、この問題は全探索か貪欲法、動的計画法でも解くことができる。

どの様に遺伝的アルゴリズムを決めるか

個体を10体とする。
適応度計算後、ソートする。
上位50%を下位50%にコピーし、コピーされた個体に交叉と突然変異を行う。
指定された世代まで行う。

今回はこれでいく。

ソースコード

精度があまり高くなく、まだ最適解以外の解も出てしまうがこれが一応遺伝的アルゴリズムを使用したコードになる

import java.util.*;

public class GA implements Algorithm{
	
	final int[] PRICE = {100,300,350,500,650};//yen
	final int[] WEIGHT = {1,2,3,4,5};//kg
	final int MAX_WEIGHT = 6;
	final int NUMBER_GENE = 8;
	final int NUMBER_ITEM = 5;
	final int MUTATE_RATE = 1;
	final char[] NAME = {'A','B','C','D','E'};
	int[][] gene;//individual
	int[][] value;//fitness and price
	
	Random rnd;
	
	static final int GENERATION = 100;
	
	public GA(){
		rnd = new Random();//init instance
	}

	@Override
	public void createIndividual() {
		gene = new int[NUMBER_GENE][NUMBER_ITEM];
		Random rnd = new Random();
		//create first generation
		for(int i = 0; i < NUMBER_GENE; i++){
			for(int j = 0; j < NUMBER_ITEM; j++){
				int init = rnd.nextInt(2);//range 0 to 1
				gene[i][j] = init;
			}
		}
	}

	@Override
	public void calcIndividual() {
		value = initValue();//initialization of this array
		int count_gene = 0,count_item = 0;
		
		for(int i = 0; i < NUMBER_GENE; i++){
			for(int j = 0; j < NUMBER_ITEM; j++){
				if(gene[i][j] == 1){
					value[count_gene][count_item] += PRICE[j];
					value[count_gene][count_item+1] += WEIGHT[j];
					
					if(value[count_gene][count_item+1] > 6){
						//if gene's weight over 6kg,this value become 0
						//and exit this loop
						value[count_gene][count_item] = 0;
						break;
					}
				}
			}
			count_gene++;
		}
	}

	@Override
	public void showIndividual() {
		for(int i = 0; i < NUMBER_GENE; i++){
			System.out.print("GENE."+(i+1));
			for(int j = 0; j < NUMBER_ITEM; j++){
				System.out.print("["+gene[i][j]+"]");
			}
			System.out.print("価値:"+value[i][0]);
			System.out.print("重量:"+value[i][1]);
			System.out.println("適応度:"+value[i][0]);// price = fitness
		}
	}

	@Override
	public void sortIndividual() {
		//value's price = fitness
		int[] sort = new int[NUMBER_GENE];
		for(int i = 0; i < NUMBER_GENE; i++){
			sort[i] = value[i][0];
		}
		
		Arrays.sort(sort);//sorted
		
		int count = 0;
		int[] array2 = new int[NUMBER_GENE];
		
		for(int i = 0; i < NUMBER_GENE; i++){
			for(int j = 0; j < sort.length; j++){
				if(sort[j] == -1)continue;
				if(value[i][0] == sort[j]){
					array2[count] = j;
					count++;
					sort[j] = -1;
					break;
				}
			}
		}
		
		int[][] result = new int[NUMBER_GENE][NUMBER_ITEM];
		int[][] result_of_value = new int[NUMBER_GENE][NUMBER_ITEM];
		
		for(int i = 0; i < NUMBER_GENE; i++){
			result_of_value[array2[i]][0] = value[i][0];
			result_of_value[array2[i]][1] = value[i][1];
			for(int j = 0; j < NUMBER_ITEM; j++){
				result[array2[i]][j] = gene[i][j];
			}
		}
		
		gene = result;
		value = result_of_value;
		
	}

	@Override
	public void selectIndividual() {
		
		int count = NUMBER_GENE/2;
		
		for(int i = 0; i < NUMBER_GENE/2; i++){
			for(int j = 0; j < NUMBER_ITEM; j++){
				gene[i][j] = gene[count][j];
			}
			count++;
		}
		
	}

	@Override
	public void crosserIndividual() {
		
		//first point and second point
		int f = 0;
		int s = NUMBER_ITEM-1;
		
		boolean flag = false;
		
		for(int i = 0; i < gene.length/2; i++){
			flag = (gene[i][NUMBER_ITEM/2] == 1)?true:false;
			if(flag){
				gene[i][NUMBER_ITEM/2] = 0;
				gene[i+gene.length/2][NUMBER_ITEM/2] = 1;
			}else{
				gene[i][NUMBER_ITEM/2] = 1;
				gene[i+gene.length/2][NUMBER_ITEM/2] = 0;
			}
		}
		
	}

	@Override
	public void mutateIndividual() {
		int random = (int)(Math.random()*10);
		for(int i = 0; i < NUMBER_GENE/2; i++){
			if(random <= MUTATE_RATE){
				if(gene[i][NUMBER_ITEM-1] == 0){
					gene[i][NUMBER_ITEM-1] = 1;
				}else{
					gene[i][NUMBER_ITEM-1] = 1;
				}
			}
		}
	}
	
	/**
	 * initialization of value -> array
	 * @return initValue<br>
	 * initialized array
	 * */
	public int[][] initValue(){
		int[][] initValue = new int[NUMBER_GENE][2];
		for(int i = 0; i < NUMBER_GENE; i++){
			for(int j = 0; j < 2; j++){
				initValue[i][j] = 0;
			}
		}
		return initValue;
	}
	
	public static void main(String... args){
		
		//make instance
		GA m = new GA();
		m.createIndividual();//create first generation
		System.out.println("第1世代");
		int count = 0;
		while(count < GENERATION){
			if(count > 0)System.out.println("第"+(count+1)+"世代");
			m.calcIndividual();
			m.sortIndividual();
			m.showIndividual();
			m.selectIndividual();
			m.crosserIndividual();
			m.mutateIndividual();
			count++;
		}
		
		System.out.println("選択したものは");
		
		for(int i = 0; i < m.NUMBER_ITEM; i++){
			if(m.gene[m.NUMBER_GENE-1][i] == 1){
				System.out.print(m.NAME[i]);
			}
		}
        System.out.println();
		
	}
}

/**
 * interface required to implements a genetic algorithm for this problem
 * */
interface Algorithm {
	
	/**
	 * method to generate the initial individual
	 * */
	public void createIndividual();
	
	
	/**
	 * method to calculate the fitness of the individual
	 * */
	public void calcIndividual();
	
	
	/**
	 * method to display the information of the generated individuals
	 * */
	public void showIndividual();
	
	
	/**
	 * method that are sorted in descending order the individual in fitness
	 * */
	public void sortIndividual();
	
	
	/**
	 * method to the selected cull the lower half of the individuals
	 * and want to copy the upper half to the lower half
	 * */
	public void selectIndividual();
	
	
	/**
	 * method to replace the random value is copied to the half in the individual
	 * */
	public void crosserIndividual();
	
	
	/**
	 * method to mutate the value ramdomly
	 * */
	public void mutateIndividual();
}

確率的解析と乱択アルゴリズム

確率的解析と乱択アルゴリズム

確率論と計算機における問題解決の手法として
確率的解析と乱択アルゴリズムについて解説していく。
ここでは以下の問題を考える。

雇用問題

問題:
あなたは秘書を雇おうとしている。
その際にあなたは代理店を利用することにした。
その代理店は毎日1人、秘書を送ってくる。
あなたはその候補者と面談し、採用するかどうか選択する。
しかし、面談するのにあなたは費用を払う必要がある。
それに加え実際に候補者を雇うには更に費用がかかる。
そこであなたは現在の秘書よりも候補者が優れている場合すぐにでも雇う。
しかし、いくら費用がかかるか見積もりたい。
かかる費用を求めなさい。

この問題における重要点

この問題の解として、候補者を逐一比較していき数値的に高ければ雇う。
といったアルゴリズムが考えられる。
またこの問題では実行時間どうこうより
計算によって算出される費用に重点が置かれるため一般的な計算量の解析とは異なりコスト解析を行う必要がある。

このアルゴリズムにおける面談コスト(C_i)は雇用コスト(C_h)より小さい。
それを踏まえてn人の候補者からm人雇用した場合
このアルゴリズムのコストはO(nC_i+mC_h)となる。
しかしn人の面談は必ず行わなければならないので、nC_iは変化しない。
よって、この問題では雇用コストC_hの解析に集中するべきである。

最悪の場合

この問題における最悪の場合は、面談した全ての候補者を雇用するときである。
つまり候補者の能力値が昇順になっているときである。
この時の雇用コストはO(nC_h)となる。

確率的解析 ~Probabilistic Analysis~

上記の問題の様に確率を用いる問題の解析を「確率的解析」という。
確率的解析とは実行時間解析の手法ではない。
上記の問題における雇用コストを解析する手法である。
確率的解析を適用するには、入力分布に関する知識を用いるか
その分布自体を仮定してしまう必要がある。
それが完了すれば実行時間の平均を計算してアルゴリズムを解析する。
平均は入力可能な分布の上でとること、つまり全ての可能な入力に対する平均を算出する。

入力分布は慎重に決める必要がある。
というのも、問題において入力集合に対してある分布をとることで効率のよいアルゴリズムを設計するためである。
しかし妥当な入力分布が存在しない場合ある。
そのようなときには確率的解析を用いることができない。

今回の雇用問題の場合、候補者はランダムに出現し二人の候補者を比較できると仮定できる。
この仮定より、以下を満たす全順序が存在する仮定する。

{ \displaystyle
集合Xにおいて関係 \le が全順序付けられるなら、Xの任意の元a,b,cに対して \\
1)a \le b かつ b \le a ならば a=b(反対称律) \\
2)a \le b かつ b \le c ならば a \le c(推移律) \\
3)a \le b または b \le aのいずれかが成立する(完全律) \\
}

ここまで示した仮定により、候補者全員を1からnまでの数でただ一通りにランク付け出来る。
それを求める関数をrank(i)とし候補iのランクを表し、大きいほど能力値が高いとする。
よって候補者がランダムに出現することと、ランクから作られるn!通りの置換のうち
任意の一つに一致する事象はどれも等確率である。
またはランクリストは一様ランダム置換と言える。

確率的解析における乱択アルゴリズム

確率的解析を用いるには入力分布に関する知識が必須となる。(上記)
しかし、入力分布がわからない場合も多く、たとえ分布がはっきりとわかっていたとしても
必ずアルゴリズムとしてモデルに出来るとは限らない。
そこでアルゴリズムの振る舞いをランダムにしてアルゴリズム設計の道具としてランダム性を使うことが出来る。
ただこのままでは最初に挙げた問題の解決法としては使えないので
問題を以下のように変更する。

代理店が候補者リストを事前に送ってくる。
そのリストから私達が毎日面談するものを決める。

こうすることでランダムな順序で送られてくる候補者に
ランダムな順序を強制する。

乱択アルゴリズム

ここまで長々と書いてきたが乱択アルゴリズムとは何か?
多少変更した上記の問題に対して、解決方がランダム性を使うものとする。
こういった場合に、アルゴリズムの振る舞いは入力データと乱数生成器の2つで決まる。
そのように乱数を使用するアルゴリズム乱択アルゴリズムという。
ただし計算機で用いる乱数生成器は、統計的にランダムに「見える」値を返す
決定アルゴリズムを用いた擬似乱数生成器によって実現されることが多い。

指標確率変数

雇用問題やその他乱択アルゴリズムを用いるであろう多数のアルゴリズム解析に指標確率変数を用いる。
これは確率と綿密な関係にある期待値をお互いに変換する便利な道具のようなものである。
標本空間Sと事象Aが与えられている時、Aに関する指標確率変数を

{\displaystyle
\begin{eqnarray}
I\{A\} = \left\{ \begin{array}{11}
1 & Aが発生する時 \\
0 & Aが発生しない時
\end{array}\right.
\end{eqnarray}
}

と定義する。

例1

一枚のコインをコイントスして表がでる回数の期待値を求めてみる。
標本空間はS=\{H,T\}とし、コインには裏表しかないので確率Pr\{H\} = Pr\{T\}=1/2となる。
今回は表が出る回数の期待値なので、事象H(表が出るとき)に対して指標確率変数X_Hを定義できる
ここでは上記の指標確率変数にしたがって、表が出るなら1、裏がでるなら0とする。

{\displaystyle
\begin{eqnarray}
X_H=I\{H\}
= \left\{ \begin{array}{11}
1 & Hが発生する時 \\
0 & Hが発生しない時(Tが発生する時)
\end{array}\right.
\end{eqnarray}
}

この様にこの問題(コイントスで表が出る回数の期待値を求める)に対して指標確率変数を定義した。
この場合の期待値をEとして、指標確率変数X_Hを用いて期待値を計算する。

{
\displaystyle
\begin{equation}
E \left[ X_H \right] \\ = E\left[ I\left\{ H \right\} \right] \\ = 1*Pr \left\{H \right\}+1*Pr \left\{T \right\} \\ = 1*(1/2)+0*(1/2) = 1/2
\end{equation}
}

よって1枚のコインをコイントスした時、表がでる回数の期待値は1/2となる。
今回のコイントス問題の様に単純な問題だと指標確率変数を用いるのは面倒に感じるかも知れないが
これが一回のコイントスではなく一般のn \in \mathbb{N}回に対してだと圧倒的に便利になる。
その例を次に示す。

例2

例1で考察した問題に対して一回のコイントスではなく一般のn \in \mathbb{N}に対して考察してみる。
ここでi番目の指標確率変数をX_iとする。この時0 \le i \le nn,i \in \mathbb{N}とする。
例1で示した様に指標確率変数の和で表せるので、全体の確率変数Xは次の様に表せる。

{\displaystyle
\begin{equation}
X = \sum_{i = 1}^{n} X_i
\end{equation}
}

上の式に対して両辺で期待値を取ると

{\displaystyle
\begin{equation}
E \left[ X \right]  \\ 
= E \left[ \sum_{i = 1}^{n} X_i \right] \\
= \sum_{i = 1}^{n} E \left[ X_i \right] \\
X_i = 1/2より \\
= \sum_{i = 1}^{n} 1/2 \\
= n/2

\end{equation}
}

となることがわかる。
これで一般の自然数nに対して期待値を指標確率変数を用いて表現することができた。
なぜ指標確率変数がすぐれた解析技術として使えるかというと、期待値が線形性を持つからである
またこの例で指した様に、指標確率変数を用いる事で大幅に計算を単純化することができる。

雇用問題への応用

候補者は上記の例同様n人いるのでX_ii人目の指標確率変数となる。
それではi人目の指標確率変数を以下のように定義する。

{ \displaystyle
\begin{eqnarray}
X_i = \left\{ \begin{array}{11}
1 & 候補者iが雇用される時 \\
0 & 候補者iが雇用されない時 \\
\end{array} \right.
\end{eqnarray}
}

また確率変数Xは次の様に定義できる。

{\displaystyle
X = X_1 + X_2 + X_3 + ... + X_n
}

更にi人目が最初の候補者からi-1人目までの候補者の中で最も優れている確率は1/iなので
上記の例2よりこの問題の期待値は次の様に表せる。

{\displaystyle
\begin{equation}
E \left[ X \right]  \\ 
= E \left[ \sum_{i = 1}^{n} X_i \right] \\
= \sum_{i = 1}^{n} E \left[ X_i \right] \\
X_i = 1/iより \\
= \sum_{i = 1}^{n} 1/i \\
= \log n +O(1)

\end{equation}
}

乱択アルゴリズムの適用

雇用問題やコイントスの問題のような入力全ての置換が等しい確率で成立するという過程が成り立つような問題では
確率的解析が乱択アルゴリズムを用いた開発の指針となる。
具体的にはアルゴリズムが実行される前に候補者の順番をランダムに並び変え
強制的に全ての置換の出現確率を等しくするようにする。

確率的解析と乱択アルゴリズムの違い

上記の記述より、候補者がランダムに出現しようと期待値は変化しない
決定的なアルゴリズム、つまり解が一意に決まるものと比べて
乱択アルゴリズムは入力時ではなく、アルゴリズム内で入力分布のランダム化が発生する。
そのため、特定の入力に対して情報の更新が何回起こるかわからない。
しかし、以前の実行と今回の実行でランダム化が同じである確率はとても低い。
よって、特定の入力があったとき、その入力に対して常に最悪の処理を引き起こすことがない
乱択アルゴリズムが最悪の処理を実行してしまうのは単に運が悪かったからというだけである。

ここまでの解説で雇用問題に乱択アルゴリズムを適用するのに必要なことは
単に入力配列のランダム化である。
これで乱択アルゴリズムが適用できる。

またこの他にも、以下のような問題にも適用できる。


追加問題

N個の箱がある。
箱を開けるまでわからないが半分は当たりになっている。
あたりを見つけよう。

この場合、最悪のケース。先頭から探索するとして
当たりが全て箱の後ろ半分に偏っていた場合O(n/2)かかってしまう。
しかし、その箱の列をランダム化すると
開ける箱の個数に対する期待値は2以下なので、すぐに発見することができる。

最後に

ここまでざっくりと乱択アルゴリズムと確率的解析について解説した
そのうち実際に乱択アルゴリズムを適用したソースコードでも載せようと思う。

gdbを使ってコアダンプの原因を解析

コアダンプは嫌いだ。

大学やその他情報系専門科のある学校に通ったことがある人が一度は触ったことがあるであろう、C言語
こいつは近年の言語に比べてものすごく面倒くさい書き方をするし、手間もかかる。
中でも最悪なのが「コアダンプ」の文字。
これは、多くのC言語ユーザを苦しめてきただろう。
一応エラーの部類に入ると思うが、こいつはどこで問題が起きているか普通は出力してくれないから直すのがすごく大変。


というわけで、今回はコアダンプの原因をgdbというデバッガを使用して解析してみる(※linuxユーザ向け)

それでは解析の準備

これはこの前自分が実際にコアダンプを発生させてしまったコード。
このファイルをここでは「quicksort.c」とする。

#include<stdio.h>

void quicksort(int array[],int left_index,int right_index){
    
    int left = left_index,right = right_index,tmp;
    
    //基点を左右の中間から選択
    int pivot = array[(left+right)/2];

    while(1){

        //左側で基点より大きい値を探す
        while(array[left] < pivot)
            left++;
        //右側で基点より小さい点を探す
        while(pivot < array[right])
            right++;//ここが間違っている

        if(right < left)
            break;

        //値の入れ替え
        //基点より小さい数値を左、大きい数値を右へ移動させる
        tmp = array[left];
        array[left] = array[right];
        array[right] = tmp;

        left++;
        right--;
    }

    //再帰関数による分割統治法を行う
    if(left_index < right)
        quicksort(array,left_index,right);
    
    if(left < right_index)
        quicksort(array,left,right_index);

}

int main(){
    int i;
    int array[] = {2,4,5,3,1,6,3,1,9};
    quicksort(array,0,9-1);
    for(i = 0; i < 9; i++)
        printf("%d\n",array[i]);

    return 0;
}

こいつを実行させると、もちろんこのメッセージがターミナルに出力される。

Segmentation fault (コアダンプ)

それでは解析を始める

まずコアダンプを解析するためにコンパイラオプションの「-g」をつける。

gcc quicksort.c -g -o quicksort
  • gオプションをつけることで、デバッグ情報を出力するようにする。

つぎに、大抵コアダンプ機能が無効にされているの以下のコマンドで有効にする

ulimit -c unlimited

それが終わったら、問題のある(コアダンプのある)実行ファイルを実行する

./quicksort
Segmentation fault (コアダンプ)

となる。

この時lsでカレントディレクトリを確認すると「core」や「core.xxxxx(xは数字)」といったファイルが生成される。


ここまでできたら、以下のコマンドでデバッグする

gdb ./quicksort core

そうすると、こんな出力が見えるはず(プログラムによって違う)

GNU gdb (Ubuntu 7.11.1-0ubuntu1~16.04) 7.11.1
Copyright (C) 2016 Free Software Foundation, Inc.
License GPLv3+: GNU GPL version 3 or later <http://gnu.org/licenses/gpl.html>
This is free software: you are free to change and redistribute it.
There is NO WARRANTY, to the extent permitted by law.  Type "show copying"
and "show warranty" for details.
This GDB was configured as "x86_64-linux-gnu".
Type "show configuration" for configuration details.
For bug reporting instructions, please see:
<http://www.gnu.org/software/gdb/bugs/>.
Find the GDB manual and other documentation resources online at:
<http://www.gnu.org/software/gdb/documentation/>.
For help, type "help".
Type "apropos word" to search for commands related to "word"...
Reading symbols from ./quicksort...done.

warning: core file may not match specified executable file.
[New LWP 13685]
Core was generated by `./base1'.
Program terminated with signal SIGSEGV, Segmentation fault.
#0  0x0000000000400616 in quicksort (array=0x7ffe0656abc0, left_index=0, 
    right_index=8) at quicksort.c:26
26	        while(pivot < array[right])
(gdb) 

これで大体わかってしまうんだけど、今回はちゃんと手順を踏む。

まずはBackTraceする。

(gdb) bt

#0  0x0000000000400616 in quicksort (array=0x7ffe0656abc0, left_index=0, 
    right_index=8) at quicksort.c:26
#1  0x0000000000400747 in main () at quicksort.c:56

これで問題のありそうな関数の番号を選択すると

#今回は#0のquicksort関数に問題がありそうだから0を選択
(gdb) frame 0

#0  0x0000000000400616 in quicksort (array=0x7ffe0656abc0, left_index=0, 
    right_index=8) at quicksort.c:26
26	        while(pivot < array[right])

というわけでこの26行目になにやら問題があることがわかった。
今回は

right++;

で配列の範囲外を参照してしまうことが問題になっていた。

というわけで、解析完了!


これがポインタなどの場合は以下の様に確かめることができる

(gdb)p [ここに確かめたい変数名、引数名など]

そうするとそれの値やメモリアドレスを指してくれる。
よくあるのは「0x0」というメモリアドレス。これはメモリが確保されていないことを指すので該当部分のポインタやアドレスまわりの処理を見なおそう。

コアダンプの原因

ここでさくっと、どうしてコアダンプが発生するかまとめておく。

  • メモリ違反
    • これが一番多いと思う
    • 配列の範囲外や未代入ポインタの参照など
  • 無限再帰や深すぎる再帰
    • 再帰から脱出できないときによくおきる

最後に

メモリ関連を特に見なおそう、配列やポインタなど。
gdbをうまく使えればコアダンプを瞬殺できるから使えると便利だけど、面倒くさいのでなるべくコアダンプが発生しないコードを書こう。

セキュリティミニキャンプ北海道2016に参加してきた

初めてミニキャンに参加してきた

LOCAL学生部に加入して、しばらくしてからミニキャンのお知らせが回ってきたから何をやるんだという感じで詳細を見てみると…。Linuxでハードウェア制御、スマホアプリでサーバサイドへの攻撃と対策、クラウドセキュリティ基礎…。これほど興味を引く内容なのに申し込まないわけにはいかない!
というわけで、申し込んで無事行くことになって、この前の週末に参加してきた。

ミニキャンに参加してみて

セキュリティあんまり興味なかったけど、セキュリティキャンプはセキュリティ以外にも色々実習を行ったり、チームに分かれてディスカッションしたりと面白かった。

印象に残った講義の感想

初日に行われた、「スマホアプリでサーバサイドアプリケーションへの攻撃と対策を学ぶ」はミニキャンプが開催された二日間の中で特に面白かった。
BurpSuiteという脆弱性診断によく使用されるソフトウェアを使用して、実際にHTTP通信を覗いたり編集したりして脆弱性を探して行ったわけだけど、チーム別に行ったのでチームのみんなで「ああだこうだ」言いながら色々試していった。別のチームの事情はよくわからないけど、自分たちのチームはセキュリティとかネットワークに強いひとがいなくてすごく苦戦した。だけど、そういう分野に強い人が居なかったが故にみんなで意見出し合いながら実習を行えたのはすごくいい経験になったし、今までセキュリティの仕事ってどんなことするのか曖昧だったものが経験を通して「これがセキュリティの仕事だ」といったものを実感できたのも良かった。


それ以外で面白かったのは「クラウドセキュリティ基礎」。
これはいろんなパワーワードが飛び交ってすごく笑えた。特に「今日のサーバはペットではなく家畜」と「耐用年数間近のサーバなんて電気代バカ食いするだけのゴミ」が最高にキテる。それにこの講義の内容はwebサービスやサーバ関連の内容からクラウドに移っていったので、最近VPSを借りてwebサービスを構築しようとしていた自分にとってはすごくためになることばかりだった。それに今までサーバとかクラウドといった環境が用意できないからという理由でweb系の事を全くやってこなかったので、本当に基礎的な部分から始まった講義は全てが新鮮だった。


その次は、なんといっても「Linuxを用いたハードウェア制御」の講義がためになった。
実はミニキャン前はこれをすごく楽しみにしていた。実際はすごくハードウェア的なことを中心にやっていくものだったけど未経験だった電子工作的な部分とシステムコール的な処理を同時に学べたのは最高に面白かった。出来ることならハードウェア制御という名前らしく、自分でハードウェアを制御するプログラムを作ってみたかった。

全体を通して

ミニキャンが開催された、二日間は最高の時間だったと振り返ってみて思う。それもそうで、普段大学では、一年生ということもあるのかこういった内容はみんな勉強しないで、たとえやったとしても「授業でやるから」「テストがあるから」といったような理由で自分がわくわくするようなことがまず無い。それに一応存在するプログラミングサークル的なものではUnityとか他のライブラリ使ってゲームを作るだけなどと、もっといろんな事を勉強して吸収していきたい自分にとってはすごく生温い環境の中にいたからである。ただ、ミニキャンは違った。ミニキャンは、みんなが何かしら自分が得られるものを全部自分のものにしてしまえ。といったような気迫ややる気と言ったものが感じられ、チームディスカッション一つとっても今までの自分にはなかった発想がたくさん展開され、演習も「そうするのか!?」といった驚きばかりだった。生温い環境に居た自分にはこういった時間がすごく刺激的だった。その中で、自分のスキルの無さを多少なりとも実感し、このままでは数年後に社会に出て、今回出会った人たちと競っていくことになったら自分は確実に劣るといった根拠の無い確証も得た。

今回、ミニキャンに参加しそういった実感があったからこそ、それをそのままにせず、自分はもっともっと成長しなければいけない、成長し続けなければいけないと思う。ミニキャンは消えそうだった自分のエンジニア魂を燃えさせるいいきっかけとなった。参加して本当に良かった。

出来ることなら来年のセキュリティキャンプの全国大会に参加したい。きっとそこにはもっと成長できる刺激的な環境が待っているだろうから。

後日談

キャンパーはその後紀伊国屋に集結してしまった。やはり似てるものはある。
あと、これからOS、カーネルの人になろうと思う。
「普段は何の人ですか?」と聞かれて答えれないのはあれだから。