有意差とは何か、シミュレーションで調べる

はじめに

いろんな分野で統計学に基づいた有意差検定が使われています。 有意差検定の方法は、さまざまな書籍、たくさんのネット記事に紹介されています。 また、RやPythonなどのツールを使えば、手軽に有意差検定をすることができます。

しかし、手持ちのデータに有意差検定を行って「有意差あり」「有意差なし」が分かったとしても、それがいったいどういうことなのか、直感的に理解するのは難しいかも知れません。

この記事では、Pythonを使ってさまざまなデータをシミュレーションで作り、さらに有意差検定を行ってみて、どういうときに有意差があり、どういうときに有意差が認められないのか、実験的に調べてみます。

Note

今回は、データの母集団が同一標準偏差正規分布に従う場合に着目します。そして、母集団から得られるサンプルの平均値に対する有意差検定を掘り下げます。

なお、本記事のPythonスクリプトの動作確認は、MacPython 2.7.15 (Anaconda)で行いました。ライブラリにnumpy, scipy, matplotlibを使っています。

例題

まず、有意差検定をする以下の例題を考えます。

問1

同じ製品を作る工場A, Bにて、それぞれ10個の製品サンプルを採取して重さを測ったところ、以下の表のとおりになった。Aの製品の重さの平均と、Bの製品の重さの平均には、有意な差があるのだろうか?

A B
99.46 g 99.29 g
100.48 g 99.17 g
103.25 g 97.63 g
98.98 g 98.14 g
99.42 g 99.14 g
100.12 g 100.56 g
100.30 g 98.73 g
100.52 g 100.12 g
100.00 g 98.94 g
101.34 g 100.33 g

問2

同じ製品を作る工場A, Bにて、それぞれ10個の製品サンプルを採取して重さを測ったところ、以下の表のとおりになった。Aの製品の重さの平均と、Bの製品の重さの平均には、有意な差があるのだろうか?

(問題文は問1と同じで、表の値が違います)

A B
99.80 g 102.01 g
100.48 g 99.70 g
99.48 g 101.27 g
99.44 g 101.23 g
101.97 g 102.35 g
101.39 g 101.89 g
100.09 g 99.00 g
100.28 g 100.63 g
100.77 g 102.67 g
101.25 g 100.56 g

回答

それでは、有意差検定をやってみます。 ここではいくつかの仮定を置きます。

  • A, Bともに、重さの母集団は正規分布に従うこととします。
  • A, Bの母集団の標準偏差は同一とします。

検定方法には、t検定のうち「独立2群の平均値の差の検定」を使います。有意水準は5%とします。

問1の回答

Pythonのscipyライブラリを使うと、手軽にt検定を行うことができます。

problem1.py

from scipy import stats
import numpy as np
# values of A
a = np.array([99.46, 100.48, 103.25, 98.98, 99.42, 100.12, 100.30, 100.52, 100.00, 101.34])
# values of B
b = np.array([99.29, 99.17, 97.63, 98.14, 99.14, 100.56, 98.73, 100.12, 98.94, 100.33])
# obtain p-value
p = stats.ttest_ind(a, b).pvalue
print "p = {}".format(p)

このスクリプトを実行すると、以下のp値が出てきます。

$ python problem1.py
p = 0.0248874728426

p値が有意水準(0.05)未満ならば、「有意差あり」と解釈されます。 よって、AとBの平均値には、有意な差があったということになります。

問2の回答

さきほどと同じスクリプトが使えますが、データは問2の値に変えます。

problem2.py

from scipy import stats
import numpy as np
# values of A
a = np.array([99.80, 100.48, 99.48, 99.44, 101.97, 101.39, 100.09, 100.28, 100.77, 101.25])
# values of B
b = np.array([102.01, 99.70, 101.27, 101.23, 102.35, 101.89, 99.00, 100.63, 102.67, 100.56])
# obtain p-value
p = stats.ttest_ind(a, b).pvalue
print "p = {}".format(p)

スクリプトを実行すると、以下のp値が出てきます。

$ python problem2.py
p = 0.181649846941

p値が有意水準(0.05)以上なので、「有意差なし」と解釈されます。 よって、AとBの平均値には、有意な差は認められなかったということになります。

実際は・・・

例題の表の値は、正規分布からサンプリングして作ったデータです。

問1の表の値は、A, Bともに、平均100, 標準偏差1の正規分布からサンプリングしたものです。 つまり、本来はAとBの母集団には差がないのに、有意差検定では「有意差あり」と判定されています。この現象は「第1種過誤」と呼ばれます。

問2の表の値は、Aは平均100、Bは平均101の正規分布からサンプリングしたものです(標準偏差はともに1)。 つまり、本来はAとBの母集団には差があるにもかかわらず、有意差検定では「有意差なし」と判定されています。この現象は「第2種過誤」と呼ばれます。

それぞれの表を作ったスクリプトも紹介しておきます。

問1の表の作成

make_problem1.py

import numpy as np
from scipy import stats

def main():
    # fix the random seed to reproduce same results 
    np.random.seed(12345)
    # set the level of significance
    sig_level = 0.05

    while True:
        # get 10 samples from a gaussian distribution (average=100, standard deviation=1),
        # and regard them as A
        a = np.random.normal(loc=100, size=(10))
        # get 10 samples from the same gaussian distribution,
        # and regard them as B
        b = np.random.normal(loc=100, size=(10))

        # do t-test between A and B
        res = stats.ttest_ind(a, b)
        # if the test result indicates "significant difference" (i.e. type 1 error), then ends.
        if res.pvalue < sig_level:
            print "type-1 error happened."
            # print sampled values of A and B
            print "       A       B"
            for i in xrange(10):
                print " %7.2f %7.2f" % (a[i], b[i])
            break

if __name__ == '__main__':
    main()

実行結果

$ python make_problem1.py
type-1 error happened.
       A       B
   99.46   99.29
  100.48   99.17
  103.25   97.63
   98.98   98.14
   99.42   99.14
  100.12  100.56
  100.30   98.73
  100.52  100.12
  100.00   98.94
  101.34  100.33
問2の表の作成

make_problem2.py

import numpy as np
from scipy import stats

def main():
    # fix the random seed to reproduce same results
    np.random.seed(12345)
    # set the level of significance
    sig_level = 0.05

    while True:
        # get 10 samples from a gaussian distribution (average=100, standard deviation=1),
        # and regard them as A
        a = np.random.normal(loc=100, size=(10))
        # get 10 samples from a gaussian distribution (average=101, standard deviation=1),
        # and regard them as B
        b = np.random.normal(loc=101, size=(10))

        # do t-test between A and B
        res = stats.ttest_ind(a, b)
        # if the test result indicates "no significant difference" (i.e. type 2 error), then ends.
        if res.pvalue >= sig_level:
            print "type-2 error happened."
            # print sampled values of A and B
            print "       A       B"
            for i in xrange(10):
                print " %7.2f %7.2f" % (a[i], b[i])
            break

if __name__ == '__main__':
    main()

実行結果

$ python make_problem2.py
type-2 error happened.
       A       B
   99.80  102.01
  100.48   99.70
   99.48  101.27
   99.44  101.23
  101.97  102.35
  101.39  101.89
  100.09   99.00
  100.28  100.63
  100.77  102.67
  101.25  100.56

どういうときに「有意差あり」「有意差なし」になるのか?

例題のように、実際には母集団に差がなくても「有意差あり」になるケース、実際には母集団に差があっても「有意差なし」になるケースがあることがわかりました。 それでは、実際にどんなときに「有意差あり」になって、どんなときに「有意差なし」になるのでしょう?

これを調べるため、次の実験をします。

実験方法

2個の正規分布A, Bから、それぞれ同個数のサンプルを採取します。

そして、Aのサンプルの平均値、Bのサンプルの平均値の間に、有意差があるかを検定します。

検定方法は、例題と同じく、t検定のなかの「独立2群の平均値の差の検定」を使います。有意水準は5%とします。

実験条件

  • Aの母集団、すなわちAのサンプルを採取する正規分布は、平均値0、標準偏差1とします。
  • Bの母集団、すなわちBのサンプルを採取する正規分布には、いろんな平均値のものを使います(0, 0.1, 1, 10)。標準偏差はつねと1にします。
  • AとBのサンプル数もいろいろと変えてみます(3〜10,000)。ただし、AとBのサンプル数はつねに同一とします。

もし、Bの母集団の平均値が0であれば、それはAの母集団と同じですので、採取されたサンプルに対しても「有意差なし」となることが期待されます。 一方、Bの母集団の平均値が0以外であれば、Aの母集団と異なりますので、「有意差あり」となることが期待されます。

母集団となる正規分布を見てみる

まず、標準偏差が1で、平均値が0〜10の正規分布がどのような形になるか、下の図で見てみましょう。

f:id:uenewsar:20181211225827p:plain

凡例"ave"が正規分布の平均値を示します。

平均0の正規分布(青)と平均10の正規分布(赤)は、山の位置が大きく異なります。そのため、これらの正規分布から採取されるサンプルの値にも大きな差が生まれますので、有意差は簡単に出そうです。

一方、平均0の正規分布(青)と平均0.1の正規分布(オレンジ)は、山の位置がほとんど重なっています。そのため、これらの正規分布から採取されるサンプルの値にも、ほとんど差がないでしょう。よって、有意差はあまり出そうにありません。

「有意差あり」になる確率を見てみる

それでは、実験をしてみます。 この実験では、AとBのサンプルの平均値に有意差が認められる確率を求めます。

それぞれの実験条件において、「Aからサンプルを採取する → Bからサンプルを採取する → 有意差検定する」というプロセスを5,000回繰り返して、そのうち何回が「有意差あり」だったかをカウントします。

そして、この「有意差あり」の回数を5,000で割ることで、AとBの平均値に有意差が認められる確率を求めます。

下の図が、その結果です。

f:id:uenewsar:20181212123341p:plain

縦軸は、有意差が認められた確率です。 横軸が、A, Bそれぞれのサンプル数です。 それぞれの線は、異なるBの母集団の平均値に対応しています。すなわち、それぞれの線で、AとBの母集団の平均値の差が異なります(凡例"diff")。

最初に、A, Bの母集団の平均値が同一の場合を見ます(青)。 期待としては、AとBの母集団は同一ですので、「有意差なし」となってほしいです。 しかし、実際には、有意水準5%の確率で「有意差あり」という間違った結果が得られることが分かります(第1種過誤)。 この第1種過誤が起こる確率は、サンプル数が多くなっても変わりません。どんなに多くのサンプルを集めようとも、設定した有意水準の確率で、間違った検定結果が得られる可能性があることを示しています。

次に、A, Bの母集団の平均値が異なる場合を見ます(オレンジ、緑、赤)。 平均値の差が10と大きい時(赤)、サンプル数によらず、ほぼ1の確率で「有意差あり」と判定されています。 しかし、平均値の差が小さくなると、サンプル数が少ない時には「有意差あり」と判定される確率が低くなっていきます。すなわち、母集団に差があるのに「有意差なし」と判定されていますので、第2種過誤が発生している状況となります。

とくに、平均値の差が0.1のときには(オレンジ)、サンプル数を相当集めないと、なかなか「有意差あり」になりません。サンプル数を10,000個集めて、ようやくほぼ1の確率で「有意差あり」という判定となります。

このグラフを作るのに使ったPythonスクリプトは以下になります。グラフが出るまで少し時間がかかります(2.3GHz Intel Core i5で2分)。

prob_of_error.py

import numpy as np
from scipy import stats
import matplotlib.pyplot as plt

def calc(num_sample, diff, sig_level):
    """
    calculate the ratio when significant differences are observed between two sets of samples.
    @param num_sample the number of samples to be picked up from a distribution.
    @param diff the difference of averages of two gaussian distributions.
    @param sig_level the level of significance.
    @return the ratio of cases when significant differences are observed between the avagages of two sets of samples.
    """

    # do 10000 trials to obtain a trustable ratio when significant differences are observed.
    num_trial = 10000
    num_significant = 0

    for i in xrange(num_trial):
        # get samples from a gaussian distribution (average=0, stdev=1),
        # and regard them as A.
        a = np.random.normal(loc=0, size=(num_sample))
        # get samples from a gaussian distribution (average=diff, stdev=1),
        # and regard them as B.
        b = np.random.normal(loc=diff, size=(num_sample))
        # do t-test between A and B.
        res = stats.ttest_ind(a, b)
        # if p-value is lower than the level of significance,
        # it means the significant difference was observed between A and B.
        # so increments the counter.
        if res.pvalue < sig_level:
            num_significant += 1

    # calculate the ratio when significant differences are observed.
    sig_ratio = float(num_significant)/num_trial

    return sig_ratio

def main():

    # fix the random seed to reproduce same results always.
    np.random.seed(12345)
    # differences of average values of distributions.
    diff_list = [0, 0.1, 1, 10]
    # numbers of samples to pick up from each distribution.
    num_sample_list = [3, 10, 30, 100, 300, 1000, 3000, 10000]

    # level of significance.
    sig_level = 0.05
    # initialize dictionary to store results.
    res = {}

    # for debug
    num_iter = len(diff_list) * len(num_sample_list)
    counter = 0

    for e_diff in diff_list:
        # initialize dictionary to store results.
        res[e_diff] = []

        for e_num_sample in num_sample_list:

            # calculate the ratio when the significant difference is observed
            # given the difference of averages of two distributions and the number of samples.
            sig_ratio = calc(e_num_sample, e_diff, sig_level)
            print "{}/{} sig_level={}, diff={}, num_sample={}, sig_ratio={}".format(counter+1, num_iter, sig_level, e_diff, e_num_sample, sig_ratio)
            
            # store the ratio.
            res[e_diff].append(sig_ratio)

            counter += 1

    # plot.
    plt.xscale("log")
    for e_diff in sorted(res.keys()):
        plt.plot(num_sample_list, res[e_diff], label="diff={}".format(e_diff), marker='o')
    plt.xlabel("number of samples [-]")
    plt.xlim((2, 30000))
    plt.ylabel("ratio when significant difference is observed [-]")
    plt.legend()
    plt.show()

if __name__ == '__main__':
    main()

実行結果

$ python prob_of_error.py
1/32 sig_level=0.05, diff=0, num_sample=3, sig_ratio=0.0501
2/32 sig_level=0.05, diff=0, num_sample=10, sig_ratio=0.0533
3/32 sig_level=0.05, diff=0, num_sample=30, sig_ratio=0.0535
        ...
31/32 sig_level=0.05, diff=10, num_sample=3000, sig_ratio=1.0
32/32 sig_level=0.05, diff=10, num_sample=10000, sig_ratio=1.0
(最後にグラフが表示される)

まとめ

この記事では、有意差がどういうときに認められ、どういうときに認められないかを直感的に理解するため、いろんなデータに対して有意差検定を行いました。今回は、データの母集団が同一標準偏差正規分布に従い、平均値のみが異なる場合に着目し、母集団から得られるサンプルの平均値に対する有意差検定を掘り下げました。

その結果、以下のことが分かりました。

  • 2群のサンプルが同一の母集団から採取されたとしても、有意水準相当の確率で、「有意差が認められる」ことがある。
  • 2群のサンプルが異なる母集団から採取されたとしても、サンプル数が少ないと、「有意差が認めらない」ケースが増える。とくに母集団の平均値の差が小さい場合、サンプル数を相当増やさないと、有意差が出ない。

タブ区切りの表をスペースを入れて整形する

エクセルなどの表計算ソフトで作った表をテキストファイルにコピーペーストしたいときがあります。

f:id:uenewsar:20151226212636p:plain

この表を全選択してコピーして、テキストエディタにペーストします。

すると、セルの間はタブ区切りになるので、表のフォーマットが崩れてしまいます。

f:id:uenewsar:20151226212638p:plain

そこで、タブ区切りの表から、表示が崩れないように整形するperlスクリプトを書きました。

タブ区切りのテキストファイル(in.txt)を標準入力で取り込むと、スペースを入れて整形して、標準出力から出力します(out.txt)。

perl makeSpaceTable.pl < in.txt > out.txt

出力されたテキストファイルは、整形されているので、表っぽく見えます。

f:id:uenewsar:20151226214100p:plain

スクリプト makeSpaceTable.pl

use strict;
use utf8;

# OSにより標準入出力の文字コードを変える
if ($^O =~ /MSWin/) {
    binmode(STDIN, ":encoding(cp932)");
    binmode(STDOUT, ":encoding(cp932)");
    binmode(STDERR, ":encoding(cp932)");
} else {
    binmode(STDIN, ":utf8");
    binmode(STDOUT, ":utf8");
    binmode(STDERR, ":utf8");
}

my @array = ();
my $maxCol = 0;

# 入力文字列をバッファに取り込む
while (<STDIN>) {
    s/[\r\n]+$//; s/\t+$//;
    if ($_ eq "") {
        next;
    }
    my @tab = split(/\t/, $_, 1000);
    push(@array, \@tab);
    # 最大カラム数を保存
    if (@tab>$maxCol) {
        $maxCol = @tab;
    }
}

# バッファの各行のカラム数を合わせる
for (my $i=0; $i<@array; ++$i) {
    for (my $j=0; $j<$maxCol-@{$array[$i]}; ++$j) {
        push(@{$array[$i]}, "");
    }
}

# 各カラムの最大文字長を求める
my @maxLen = ();
for (my $i=0; $i<$maxCol; ++$i) {
    my $tmp = 0;
    for (my $j=0; $j<@array; ++$j) {
        if (getLenForPrint($array[$j][$i]) > $tmp) {
            $tmp = getLenForPrint($array[$j][$i]);
        }
    }
    push(@maxLen, $tmp);
}

# 表示する
for (my $i=0; $i<@array; ++$i) {
    my $buf = "";
    for (my $j=0; $j<$maxCol; ++$j) {
        # カラム間に挿入するスペースを作る
        my $space = " " x ($maxLen[$j] - getLenForPrint($array[$i][$j]));
        if ($i==0 || $j==0) {
            # 第1行、または第1カラムの場合、左詰め
            $buf .= $array[$i][$j] . $space;
        } else {
            # それ以外は右詰め
            $buf .= $space . $array[$i][$j];
        }
        # カラム間にはスペース3個を入れる
        $buf .= "   ";

    }
    # 行末の空白は消しておく
    $buf =~ s/ +$//;
    # 出力
    print $buf . "\n";
}


# 表示される文字長さを数える
# (半角=1, 全角=2になるように)
sub getLenForPrint {

    my ($in) = @_;

    # 半角文字だけを残す
    my $ascii = $in;
    $ascii =~ s/[^\p{InBasicLatin}]//g;

    # 非半角文字だけを残す(→全角文字とみなす)
    my $nonAscii = $in;
    $nonAscii =~ s/[\p{InBasicLatin}]//g;

    return (length($ascii) + length($nonAscii) * 2);

}

ファミコンのDPCMで再生した音をシミュレートする

ファミコンの音源には、サンプリング音声を再生するDPCMチャンネルがあります。 用意した音声ファイルをサンプリングしてDPCMチャンネルから再生したときに、どのような音になるかをシミュレートするpythonスクリプトを書きました。

DPCMチャンネルの仕様は、こちらの記事を参考にしました。
FC音源とは (ファミコンオンゲンとは) [単語記事] - ニコニコ大百科

動作確認は、Macpython 2.7.6 (homebrew)で行いました。ライブラリにnumpyを使っています。

スクリプト

dpcm.py

# coding: UTF-8

# 入力の音声ファイルから、
# あたかもRP2A03のDPCMで再生したような音声ファイルを作る

import wave
import sys
import struct
import numpy as np
import math

# 音声データをwavファイルとして保存する
# 音声データは、モノラル, 量子化ビット数16とする。
# @param inData 音声データ
# @param inSampFreq サンプリング周波数[Hz]
# @param inFileName 出力するファイルのファイル名
def writeWav(inData, inSampFreq, inFileName):
    data = struct.pack("h" * len(inData), *inData)
    wf = wave.open(inFileName, "w")
    wf.setnchannels(1)
    wf.setsampwidth(2)
    wf.setframerate(inSampFreq)
    wf.writeframes(data)
    wf.close()

# 音声ファイルを読み込む
# 非圧縮, 1サンプル2バイト(量子化ビット数16)のwavファイルを読み込む。
# チャンネル数は1または2。チャンネル数が2の場合には、Lチャンネルのみを読み込む。
# @param inFileName 音声ファイル名(.wav)
# @return [音声データ, サンプリング周波数[Hz]] 
def readWav(inFileName):

    # 音声ファイルを開く
    wr = wave.open(inFileName, "rb")

    # 音声ファイルの検査
    # 1サンプルのバイト数は2であること
    if wr.getsampwidth() != 2:
        print " error : size per one sample must be 2."
        wr.close()
        usage()
        quit()
    # 非圧縮であること
    if wr.getcomptype() != "NONE":
        print " error : input file must be uncompresed format."
        wr.close()
        usage()
        quit()
    # チャンネル数は1か2であること
    channels = wr.getnchannels()
    if channels != 1  and channels != 2 :
        print " error : channel number must be 1 or 2."
        wr.close()
        usage()
        quit()

    # 読み込む
    sampFreq = wr.getframerate()
    data = wr.readframes(wr.getnframes())
    numData = np.frombuffer(data, dtype=np.int16)
    wr.close()

    # 読み込んだデータをintのリストに付け替える
    i=0
    out = []
    while i<len(numData):
        out.append(int(np.asscalar(np.float64(numData[i]))))
        # 音声ファイルが2チャンネルならLだけを取る
        i += channels

    return [out, sampFreq]

# DPCMチャンネルで再生した音を作る
# 参考URL: http://dic.nicovideo.jp/a/fc音源
# ロジック:
# 1. 圧縮後サンプル値を作る
#   1個目のサンプル
#     元のサンプル値(16 bit)を、7 bitに圧縮する
#   2サンプル以降:
#     当該サンプル値(16 bit)を、7 bitに圧縮する
#     直前の圧縮後サンプル値(7 bit)と比較し、
#     ・当該サンプルのほうが大きければ、直前の圧縮後サンプル値+1を追加する
#     ・当該サンプルのほうが小さいまたは等しければ、直前の圧縮後サンプル値-1を追加する
# 2. 音声に変換する
#   サンプル値の下位1 bitを0にする
#   (デルタ初期ボリュームは7bit長だが、再生時は下位1bitを無視した64段階で
#    ボリューム変更される)
#   そのあと、9 bit(=16-9)だけ左シフトし、量子化ビット数を16に戻す。
# @param inData 入力となる音声波形
# @param inSampFreq 入力音声のサンプリング周波数 [Hz]
# @param inRegId DPCM化後のサンプリング周波数を決める値
def convDpcm(inData, inSampFreq, inRegId):

    # DPCM再生のサンプリング周波数を決める
    # 指定された0-15の値に対応する値をテーブルから引く
    # サンプリング周波数[Hz] = 1789772.5 / テーブルから引いた値
    if inRegId < 0 or inRegId > 15:
        print " error : number for determining sampling frequency must be from 0 to 15."
        usage()
        quit()
    dpcmFreqTable = [428,380,340,320,286,254,226,214,190,160,142,128,106,85,72,54]
    dpcmFreq = 1789772.5 / dpcmFreqTable[inRegId]
    print " DPCM-converting sampling frequency[Hz]:", dpcmFreq

    # 音声の総時間
    totalTime = float(len(inData)) / inSampFreq
    
    # 初期値の設定
    # 注目サンプルの時刻
    t = 0.0
    # 直前サンプルの圧縮後の値
    prev = 0.0
    # 出力
    out = []

    while t < totalTime:

        # 注目サンプルの値を取る
        sample = inData[int(math.floor(t * inSampFreq))]

        if t <= 0.0:
            # 第1サンプルのとき
            # 音を7 bitに圧縮する
            # pythonは算術シフトなので、普通に9 bit(=16-7)だけ右シフトさせればよい
            val = sample >> 9
        else:
            # 第2サンプル以降
            # 当該サンプルを圧縮する
            tmpVal = sample >> 9
            # 圧縮後の当該サンプルが直前のサンプルより大きければ、直前の圧縮後サンプルの+1にする。
            # 以下ならば-1する。
            if tmpVal > prev:
                val = prev+1
            else:
                val = prev-1

        # 圧縮後サンプル値が範囲を出ていれば、範囲内に収める
        # 範囲は、-64〜63 (7 bitの符号付き整数)
        if val > 63:
            val = 63
        elif val < -64:
            val = -64

        # 出力のリストに登録
        out.append(val)
        # 当該圧縮後サンプル値は次につかうので、保存
        prev = val

        # 時刻を更新
        t += 1.0 / dpcmFreq

    # 圧縮後サンプル値を、量子化ビット16 bitの波形に戻す
    # 圧縮後サンプル値は7 bitだが、下位1bitはボリュームは再生時のボリュームには反映されない(無視される)ため、
    # 下位1 bitを削除した後に、16 bitに戻す
    for i in xrange(len(out)):
        out[i] = (out[i] >> 1) << 10

    # 出力を返す
    return [out, dpcmFreq]

# usage
def usage():

    print "usage: python dpcm.py (input file) (output file) [options]"
    print " This script converts input speech so as to be played from NES DPCM channel."
    print "  (input file) : input wav file name."
    print "  (output file) : output wav file name."
    print "  [options]"
    print "   -ds (number) : specify DPCM sampling frequency."
    print "     The specified number must be in the range from 0 to 15 [defalut=15]."
    print "     Each number indicates a stored value to a register that determines"
    print "     DPCM sampling frequency."
    print "     The relationship between the specified numbers and the sampling"
    print "     frequencies are as followings:"
    print "       (number) (frequency[Hz])"
    print "             0            4182"
    print "             1            4710"
    print "             2            5264"
    print "             3            5593"
    print "             4            6258"
    print "             5            7046"
    print "             6            7919"
    print "             7            8363"
    print "             8            9420"
    print "             9           11186"
    print "            10           12604"
    print "            11           13983"
    print "            12           16885"
    print "            13           21056"
    print "            14           24858"
    print "            15           33144"


# メインルーチン
if __name__ == "__main__" :

    # 引数を取る
    if len(sys.argv) < 3:
        print " error : less args."
        usage()
        quit()
    
    # 入力ファイル名を受け取る
    inFileName = sys.argv[1]
    # 出力ファイル名を受け取る
    outFileName = sys.argv[2]

    # DPCMのサンプリングレートを決める値を受け取る
    # (デフォルト値)
    dpcmFlag = 15
    i=3
    while i<len(sys.argv):
        if sys.argv[i] == '-ds':
            i+=1
            if i>=len(sys.argv):
                print " error : invalid args."
                usage()
                quit()
            if not sys.argv[i].isdigit():
                print " error : invalid args."
                usage()
                quit()
            dpcmFlag = int(sys.argv[i])
        else:
            print " error : invalid args."
            usage()
            quit()
        i+=1
        
    # 入力音声ファイルを受け取る
    [data, sampFreq] = readWav(inFileName)

    print " input file:", inFileName
    print " output file:", outFileName
    print " input sampling frequency [Hz]:", sampFreq
    print " speech length [s]:", float(len(data)) / sampFreq

    
    # DPCM変換 → 復元を実行
    [out, dpcmFreq] = convDpcm(data, sampFreq, dpcmFlag)

    # 出力時のサンプリング周波数
    # wavの場合、サンプリング周波数は整数になるので、
    # DPCM変換時の周波数の四捨五入にする
    outSampFreq = round(dpcmFreq)
    print " output sampling frequency[Hz]:", outSampFreq

    # できた音声ファイルを出力する
    writeWav(out, outSampFreq, outFileName)

使い方

用意した音声ファイル(in.wav)から、DPCMチャンネルからの再生をシミュレートした音声ファイル(out.wav)を作るには、以下のコマンドを実行します。

$ python dpcm.py in.wav out.wav

ファミコンDPCMチャンネルでは、再生音のサンプリング周波数を16個のなかから選択できます。再生音のサンプリング周波数は、デフォルト設定でもっとも高い33144Hzになります。再生音のサンプリング周波数を変更するには、-dsオプションで値(0-15)を指定します。この値とサンプリング周波数の対応は、スクリプトのusageメッセージに記載しています。

実際のファミコンでは、再生音のサンプリング周波数は小数になるようです。このスクリプトでは、wavファイルとして出力する関係上、再生時のサンプリング周波数は、整数に四捨五入しています。

例として、ドラム音を使って、DPCMチャンネルから再生した音を作ります。再生音のサンプリング周波数の設定は、デフォルト(33144Hz)です。

元のドラム音 Loop01 by uenewsar on SoundCloud - Hear the world’s sounds

DPCMチャンネルからの再生音に変換した音 Converted by uenewsar on SoundCloud - Hear the world’s sounds

ドラム音には、以下のサイトのものを使わせていただきました。
Music is VFR ~BGM・ジングル・効果音のフリー素材~