量子情報理論の基本：データ圧縮（１）

$$
\def\bra#1{\mathinner{\left\langle{#1}\right|}}
\def\ket#1{\mathinner{\left|{#1}\right\rangle}}
\def\braket#1#2{\mathinner{\left\langle{#1}\middle|#2\right\rangle}}
$$

## はじめに

[前々回の記事](https://qiita.com/SamN/items/921c529c0eb09867c5cb)で「ホレボー限界」を取り上げ、量子通信チャネルの限界について説明しました。そのときの想定は、古典情報を量子ビットで符号化して送信し、受信側でPOVM測定することで古典情報を復号するという、「古典-量子-古典」という通信路のお話でした。今回、取り上げたいのは、そうではなくて、量子状態そのものを「情報」として伝送するお話です。古典的なシャノンの情報理論において、情報源の性質（確率分布）がわかったら、最低限必要な通信容量が決まるという有名な定理がありますが、その定理の量子版がシューマッハによって証明されています。そのあたりの理解まで何とか行きたいと思っています。ちょっと長くなりそうなので、２回に分けます。１回目の今回は、まず前提知識として、古典情報理論における「雑音のないチャネルにおけるシャノンの符号化定理」について説明し、２回目の次回に、本題である「雑音のない量子チャネルにおけるシューマッハの符号化定理」を説明します。各回、理解ができたところで、Pythonプログラムで、その符号化動作をシミュレーションしてみます[^1]。

[^1]: ちなみに、いま「符号化」という言葉を使いましたが、「データ圧縮」と同義と思っていただいて良いです。「符号化」というのは情報をビット列に変換することです。その際、できるだけ効率的に、つまり必要なリソースを少なくして「通信」したり「蓄積」したいわけです。そのため、情報をできるだけ「圧縮」したいというモチベーションが生まれ、シャノン以降、長年研究がなされてきました。ちょっとうるさいことを言うと、正確には、「データ圧縮」は「符号化」の目的の一つだと思うのですが、非常に大きな目的なので、色々な文脈の中で、ほぼ同じ意味で使われているというのが現状かと思います。また、今回、「雑音のない」という限定をつけていますが、当然、「雑音のある」場合もあります。未勉強なので、理解できたら機会を改めて、紹介したいと思います（予定は未定です）。

参考にさせていただいたのは、以下の文献です。

1. [ニールセン、チャン「量子コンピュータと量子通信（３）」オーム社（2005年）](https://www.ohmsha.co.jp/book/9784274200090/)
2. [石坂、小川、河内、木村、林「量子情報科学入門」共立出版（2012年）](https://www.kyoritsu-pub.co.jp/bookdetail/9784320122994)
3. [富田「量子情報工学」森北出版（2017年）](https://www.morikita.co.jp/books/book/3109)

## i.i.d.情報源

これから考えていきたいのは、一言で言うと、情報源の性質がこうだったときに、原理的にどれだけのデータ圧縮ができるだろうかという疑問への答えです。そのため、情報というある種つかみどころのないものを抽象化して、モデル化して、数学的に扱える形にしておく必要があります。本質を捉えるための単純化によって、現実の情報から多少外れてしまう可能性もありますが、実用上有益な結果が得られるのであれば、それも良しとします。という考え方で、まず、情報源というものを抽象化・モデル化します。

例えば、アルファベット文字列からなる英文のテキスト情報をイメージしてみてください。"I have a dream."という英文があったときに、アルファベット文字は"I","$\space$","h","a","v","e"...という具合に、時系列順に生起するという見方ができます。世の中に英文というのは、これ以外にも大量にあって、それらを見ると生起しやすい文字（"e"とか"a"とか）とあまり生起しない文字（"q"とか"z"とか）がありますので[^2]、各時点で生起するアルファベット文字を確率変数と考えてみるというのが、一つの抽象化・モデル化のアイデアになります[^3]。このアルファベット文字に関する確率は、どんな時点でも共通であるというのは自然な仮定としてあり得そうですし、さらに、異なる時点で発生したアルファベット文字同士は、互いに無関係に発生するという仮定を置いても、良さそうです[^4]。このように、情報源から別々に発生した値は互いに独立(independent)で、同一分布(identically distributed)に従うという仮定をおいた情報源のことを「i.i.d.情報源(independent and identically distributed information source)」と呼びます。これを前提に以降の議論を進めていきます。

[^2]: 参考：[文字頻度表](http://www7.plala.or.jp/dvorakjp/hinshutu.htm)

[^3]: "I have a dream."というのは、マーチン・ルーサー・キング・ジュニアの有名なスピーチの一節で、英語の授業で読んだり、聞いたりした人は多いと思います。何度聞いても、そのスピーチは感動的で、彼が何を言わんとしていたのか、この"I have a dream."という一節を見るだけで、伝わってくるものがあります。そういった一切の含意を捨てて、抽象化して、ひたすら数学的な対象として、情報を捉えましょう、ということをこれから説明していくわけです。シャノンがはじめた情報理論がそれです。とてもクールな立場ですが、これによって世界は大きく変わりました。スピーチの動画、リンク張っておきます。"I have a dream."を連呼し始めるのは、11分40秒あたりからです。[【日本語字幕】キング牧師演説"私には夢がある" - Martin Luther King "I Have A Dream"](https://www.youtube.com/watch?v=eQ6q2cnVXqQ)

[^4]: 英文をイメージすると、正確には違います。例えば、"t"の後に”h”が来る確率は、他の場合と比べて多いです。が、多くの情報源おいて、大雑把にこんな仮定をしても、そう外れた仮定ではなさそうです。また、もっともシンプルな仮定からはじめて、次第に複雑なものを考えていくのが、理論構築をする常套手段です。というわけで、まずこのようなシンプルな仮定から議論をはじめます。


## 典型系列

シャノンの符号化定理の前に、「典型系列」「非典型系列」という概念について説明しておく必要があります。いま、i.i.d.情報源が確率変数$X_1,X_2,X_3, \cdots$を作り出すとします。例えば、コインを投げて表が出るか裏が出るかという事象を確率変数と見立て、それを何度も繰り返して表裏の系列つくるようなことをイメージしてもらえば良いです。表の場合"0"、裏の場合"1"という数値を割り当てると、その系列は２進数の系列になります。普通のコインであれば表も裏もどちらも確率1/2で生起しますが、なるべく一般的に議論したいので、表が出る確率が$p$であるように作られているコインであるとしておきます。という想定だとしても、毎回同じコインを投げるので、「同一分布(identically distributed)」であるというのは自然な仮定です。また、どっちが出るかは、過去の結果に依存しない（例えば、前回「表」だったので今回「表」が出やすくなる、わけではない）ので、「独立(independent)」というのも変な仮定ではありません。というわけで、コイン投げで得られる情報は、i.i.d.情報源と見なしても良いです。

このi.i.d.情報源の無限の系列の中から$n$個の確率変数$X_1,X_2,\cdots,X_n$を取り出してみると、その実際の値$x_1,x_2,\cdots,x_n$にはいろんな可能性があるのですが、生起しやすい系列パターンとめったに生起しない系列パターンに分けることができます。生起しやすい系列を「典型系列(typical sequence)」、めったに生起しない系列を「非典型系列(atypical sequence)」と呼びます。0が出る確率を$p=1/4$として、取り出す系列の数を$n=100$とします。そうすると、0が出る回数は、だいたい25回あたりになることが想像できます。この100個のサンプリングを何度も繰り返したとしても、おそらく多少の揺れはあるでしょうか、0が出る回数は、25回あたりに集中するはずです。逆に、5回しかないとか、85回もあったら、これは非常に珍しい現象が起きたか、コイン投げをする人がイカサマ師であるかのどちらかです。いまの例で言うと、0が出る回数が25回くらいになっている系列のことを「典型系列」と言います。

典型系列は生起しやすい系列であると言いいましたが、では、その確率はどのように定量化できるでしょうか。確率変数系列の中から、n個のサンプリングをして、系列$x_1,x_2,\cdots,x_n$を得る確率を$p(x_1,x_2,\cdots,x_n)$とすると、i.i.d.情報源なので、

```math
p(x_1,x_2,\cdots,x_n) = p(x_1)p(x_2) \cdots p(x_n) \approx p^{np} (1-p)^{n(1-p)}  \tag{1}
```

となります。両辺の対数をとると、

```math
-\log p(x_1,x_2,\cdots,x_n) \approx -np \log p - n(1-p) \log (1-p) = n H(X)  \tag{2}
```

なので、

```math
p(x_1,x_2,\cdots,x_n) \approx 2^{-nH(X)}  \tag{3}
```

という具合にエントロピーを使って近似できます。ここでは、確率変数が"0"と"1"という値しかとらない２値をイメージして、式(3)を導きましたが、多値の場合でも成り立ちます。証明してみます。

【証明】

確率変数$X_i \space (i=1,2,\cdots)$において生起するシンボルが$A=\\{ a_1,a_2,\cdots,a_{|A|}\\}$のいずれかだったとします。$|A|$は集合$A$の要素数を表します。確率変数系列からn個をサンプリングし、その値が$\\{ x_1,x_2,\cdots,x_n\\}$だった確率は、i.i.d.情報源を前提としているので、以下のように、各確率の積で表すことができます。

```math
p(x_1,x_2,\cdots,x_n) = p(x_1)p(x_2) \cdots p(x_n)  \tag{4}
```

また、$p(x_i)$の関数形は同一です。$x_i=a_{\mu}$のときの確率$p(a_{\mu})$は、$n$個の系列の中で$a_{\mu}$が登場する回数を$N(a_{\mu})$と書くとすると、

```math
p(a_{\mu}) \approx \frac{N(a_{\mu})}{n}  \tag{5}
```

と近似できます。そうすると、

```math
p(x_1,x_2,\cdots,x_n) \approx \prod_{\mu=1}^{|A|} p(a_{\mu})^{N(a_{\mu})} = \prod_{\mu=1}^{|A|} p(a_{\mu})^{np(a_{\mu})}  \tag{6}
```

となります。両辺の対数をとると、

```math
- \log p(x_1,x_2,\cdots,x_n) \approx \sum_{\mu=1}^{|A|} np(a_\mu) \log p(a_{\mu}) = nH(X) \tag{7}
```

なので、

```math
p(x_1,x_2,\cdots,x_n) \approx 2^{-nH(X)}  \tag{8}
```

となり、式(3)が同様に証明できました。（証明終）

確率の逆数が、それが生起する場合の数なので、典型系列の数はだいたい$2^{nH(X)}$くらいと思っておけば良いです。ということは、情報源をn個のブロックに区切って符号化することを考えると、$nH(X)$ビットの符号を用意しておけばほとんどの情報は間違いなく伝送できるということになります。もちろん、非典型系列が発生する場合もあるので、その場合は、諦めて適当な符号を割り当てます。とすると、一定の割合で伝送誤差が発生しますが、nを十分に大きくすると、その割合を十分に小さくすることができて、結局、信頼できる通信を実現することが可能になります。シャノンが証明したのは、まさにこのことです。この証明に入る前に、典型系列について、もう少し、深い理解が必要なので、それについて説明します。

$\epsilon > 0$が与えられたとき、

```math
2^{-n(H(X)+\epsilon)} \leq p(x_1,x_2,\cdots,x_n) \leq 2^{-n(H(X)-\epsilon)}  \tag{9}
```

を満たす系列$x_1,x_2,\cdots,x_n$を「$\epsilon$典型系列」と言います。長さ$n$のそのようなすべての$\epsilon$典型系列の集合を$T(n,\epsilon)$と表します。式(9)は、

```math
| \frac{1}{n} \log \frac{1}{p(x_1,x_2,\cdots,x_n)} - H(X) | \leq \epsilon  \tag{10}
```

のように書きかえられます。これで、漠然と定義してきた典型系列を数学的に扱えるように一応定量化できました。ということで、次に、典型系列に関する３つの定理について順に説明していきます。

### 典型系列の定理（１）

[ニールセン、チャン](https://www.ohmsha.co.jp/book/9784274200090/)の記載をそのまま引用します。

「$\epsilon > 0$を固定する。任意の$\delta > 0$、十分大きな$n$に対して文字列が$\epsilon$典型的である確率は少なくとも$1-\delta$である。」

さて、何を言っているのかさっぱりわかりません、という人も多いかと思うので、少し噛み砕きます。$\epsilon$というのは、典型系列の範囲を決める$\epsilon$です。いまこの$\epsilon$を固定したとします。頭の中で何か適当な正の数字（例えば、0.1とか）を思い浮かべてください。で、もうひとつ$\delta$というのが登場します。これは任意の正の値として設定しても良いです。この定理が主張しているのは、どんな$\delta>0$を決めたとしても、十分に大きな$n$をとれば、長さ$n$の系列が、$\epsilon$典型系列になる確率は、少なくとも$1-\delta$である、ということです。例えば、$\delta$を割と大きい値（0.999とか）にすると、この定理の主張はかなり自明な感じがします。なぜなら、$\epsilon$典型になる確率はとても小さい＝とてもゆるい条件になるので、そんなに大きな$n$でなくても成り立ちそうな気がします。一方、$\delta$の値が非常に小さい値だとするとどうでしょうか（例えば、0.000001とか）？この定理の主張は非常に厳しい条件になります。なぜなら、$\epsilon$典型になる確率は、ほぼ1であると言っているので。しかし、$n$の値は十分に大きな値をとれます。十分大きな$n$をとれば、$\epsilon$典型になる確率を$1-\delta$以上にすること（すなわち十分1に近づけること）ができます。ということを、この定理は言っています。ざっくり一言で言うと、とても簡単な話で、十分大きな$n$をとれば、どんな系列も$\epsilon$典型系列になるということです。サイコロを転がして出る目を記録していく実験を考えてみてください。10回転がした程度では、典型系列になっていなかったとしても、100回とか1000回とか転がせば、きっと典型系列になっているに違いありません。ということを数学的に厳密に言っているだけですね。では、証明してみます。

【証明】

$X_1,X_2,\cdots$はi.i.d.情報源から得られる確率変数系列とします。$p(X_i)$は独立で同一分布なので、$-\log p(X_i)$も独立で同一分布です。ということは、$-\log p(X)$の期待値$E(-log p(X))$は、十分大きな$n$個の$-\log p(X_i)$を取り出したときの平均値に近似することができます。すなわち、任意の$\epsilon > 0, \delta > 0$に対して、

```math
p(| \sum_{i=1}^{n} \frac{- \log p(X_i)}{n} - E(- \log p(X)) | \leq \epsilon) \geq 1-\delta  \tag{11}
```

が成り立つ$n>0$が必ず存在します。これは大丈夫でしょうか。どんなに小さな$\epsilon$や$\delta$をとったとしても、十分に大きな$n$をとれば、上の式を満たすことができます、ということです。つまり、十分大きな$n$をとれば$E(- \log p(X)) = \sum_{i=1}^{n} \frac{- \log p(X_i)}{n}$にできるということで、これは、いわゆる「大数の法則」の言い換えに過ぎません。

ここで、

```math
\begin{align}
& E(- \log p(X)) = H(X) \\
& \sum_{i=1}^{n} \log p(X_i) = \log(p(X_1,X_2,\cdots,X_n))  \tag{12} 
\end{align}
```

なので、式(11)は、

```math
p(| - \frac{\log(p(X_1,X_2,\cdots,X_n))}{n} - H(X) | \leq \epsilon) \geq 1-\delta  \tag{13}
```

となります。この左辺の$p$の中身は$\epsilon$典型系列の定義そのものなので、式(13)は、$\epsilon$典型系列である確率が少なくとも$1-\delta$であることを言い表しています。ということで、証明ができました。（証明終）


### 典型系列の定理（２）

[ニールセン、チャン](https://www.ohmsha.co.jp/book/9784274200090/)の記載を引用します。

「任意の固定した$\epsilon >0$、$\delta >0$、十分大きな$n$に対して$\epsilon$典型系列の数$|T(n,\epsilon)|$は次式を満たす。

```math
(1-\delta) 2^{n(H(X)-\epsilon)} \leq |T(n,\epsilon)| \leq 2^{n(H(X)+\epsilon)}  \tag{14}
```
」

ざっくり言い換えると、$\epsilon \to 0$、$\delta \to 0$の極限で、十分大きな$n$をとれば、$|T(n,\epsilon)| \approx 2^{n(H(X)}$に近似できる、ということです。先程、式(8)の下あたりで「典型系列の数はだいたい$2^{nH(X)}$くらい」と言いましたが、そのことを数学的に厳密に表現すると、こんな定理になります。では、証明します。

【証明】

系列$x=(x_1,x_2,\cdots,x_n)$が生起する確率に対する式(9)と確率が1を超えないことから、

```math
1 \geq  \sum_{x \in T(n,\epsilon)} p(x) \geq \sum_{x \in T(n,\epsilon)} 2^{-n(H(X)+\epsilon)} = |T(n,\epsilon)| 2^{-n(H(X)+\epsilon)}  \tag{15}
```

なので、

```math
|T(n,\epsilon)| \leq 2^{n(H(X)+\epsilon)}  \tag{16}
```

が成り立ちます。これで、式(14)の右側の不等号が証明できました。次に左側の不等号についてです。

典型系列の定理（１）および式(9)より、

```math
1-\delta \leq \sum_{x \in T(n,\epsilon)} p(x) \leq \sum_{x \in T(n,\epsilon)} 2^{-n(H(X)-\epsilon)} = |T(n,\epsilon)| 2^{-n(H(X)-\epsilon)}  \tag{17}
```

なので、

```math
(1-\delta) 2^{n(H(X)-\epsilon)} \leq |T(n,\epsilon)|  \tag{18}
```

が成り立ちます。式(16)と式(18)を合わせると、式(14)になります。（証明終）


### 典型系列の定理（３）

[ニールセン、チャン](https://www.ohmsha.co.jp/book/9784274200090/)の記載を引用します。

「$S(n)$は情報源からの長さ$n$の文字列の集まりで、その大きさはたかだか$2^{nR}$であるとする。ここで$R\<H(X)$は固定されている。このとき任意の$\delta\>0$と十分大きな$n$に対して次式が成り立つ。

```math
\sum_{x \in S(n)} p(x) \leq \delta  \tag{19}
```

」

典型系列の数は、近似的に$2^{nH(X)}$くらいですが、$H(X)$よりも小さい$R$を使って、$2^{nR}$個くらいしかない部分集合$S(n)$を考えるわけです。$n$を十分に大きくしていくと、長さ$n$の系列集合の中に含まれる$S(n)$の割合が十分に小さくなるということを主張しています。$n \to \infty$の極限で、その割合は0に収束します。ホントですか？と言いたくなるかもしれませんが、ホントです。以下で証明してみます。2のべき乗の威力をご堪能ください。

【証明】

$S(n)$の文字列を典型系列と非典型系列に分けて考えます。典型系列の定理（１）より、$n$が十分に大きいとき、長さ$n$の系列集合に含まれる非典型系列の数は十分に小さくなるので、その部分集合である$S(n)$に含まれる非典型系列の数も十分に小さくなります。$n \to \infty$の極限では0になります。

一方、$S(n)$の中に含まれる典型系列の数は、$S(n)$に含まれる系列数以下です。つまり、たただか$2^{nR}$です。典型系列の出現確率は、$2^{-nH(X)}$だったので、$S(n)$が出現する確率は、たかだか、$2^{n(R-H(X))}$となります。つまり、十分大きな$n$に対して、

```math
\sum_{x \in S(n)} p(x) \leq 2^{n(R-H(X))}  \tag{20}
```

です。いま、$R < H(X)$だったので、任意の$\delta > 0$に対して


```math
\sum_{x \in S(n)} p(x) \leq \delta  \tag{21}
```

を満たす、$n>0$が存在するということになります。つまり、$n \to \infty$の極限で、

```math
\sum_{x \in S(n)} p(x) \to 0  \tag{22}
```

が成り立ちます。（証明終）


## 雑音のないチャネルにおけるシャノンの符号化定理

さて、それでは本日のメイン・イベントであるシャノンの定理を証明する準備が整いました。またまた、[ニールセン、チャン](https://www.ohmsha.co.jp/book/9784274200090/)の記載を引用します。

「$\\{ X_i \\}$はエントロピー$H(X)$のi.i.d.情報源であるとする。$R>H(X)$とする。このときこの情報源に対してレート$R$の信頼できる圧縮法が存在する。逆にもし$R<H(X)$ならば如何なる圧縮法も信頼性が高くない。」

一言で言うと、情報源のエントロピーよりも、ちょっとでも大きい伝送レートを用意しさえすれば、この情報を正確に伝送する圧縮法を構築することができるということです。では、証明してみます。

【証明】

まず、$R>H(X)$の場合です。

$H(X)+\epsilon \leq R$を満たすような$\epsilon$を選び、その$\epsilon$で$\epsilon$典型系列の集合$T(n,\epsilon)$を考えます。典型系列の定理（１）と（２）および$R>H(X)$という前提より、任意の$\delta>0$と十分大きな$n$に対して、

```math
|T(n,\epsilon)| \leq 2^{n(H(X)+\epsilon)} \leq 2^{nR}  \tag{23}
```

が成り立ちます。つまり、典型系列の個数は、たかだか$2^{nR}$個で、そのような系列を生起させる確率は、少なくとも$1-\delta$です。もう少し噛み砕いて言うと、十分に大きな$n$をとれば、その系列は確率$1$で典型系列となり、かつ、その個数は$2^{nR}$個です[^5]。従って、$n$個からなる系列を$nR$ビット使って符号化すれば、十分正確に伝送できるということになります。これを称して「レートRの信頼できる圧縮法が存在する」という言い方をしています。

[^5]: $\epsilon$とか$\delta$が出てくる数学の言い回しは独特で最初面食らうと思いますが、厳密な言い方をしようとすると、どうしてもこうなってしまいます。が、よく吟味すると、意外とやさしいことを言っていることも多々あります。今の場合、サイコロを十分に何度も転がせば、それは典型系列になっていて、かつ、その典型系列の個数はたかだか$2^{nR}$個（ただし$R$はサイコロ１回転がしの場合のエントロピーよりも少し大きい値）になりますと言っています。

次に、$R<H(X)$の場合です。

典型系列の定理（３）より、$R<H(X)$のとき、情報源から出力された系列が$2^{nR}$個の系列集合$S(n)$に含まれる確率は、十分大きな$n$に対して$0$になります。従って、正確に伝送できる圧縮方法は存在しません。（証明終）


## 符号化シミュレーション

情報源のエントロピーに相当する伝送レート近くにまで可逆にデータを圧縮する方法として、「ハフマン符号」とか「算術符号」とか色々提案されていますが、ここでは、シャノンの定理を実感してみることのみを目的として、非常に原始的な符号化を実装することにします。すなわち、十分に大きな$n$をとれば、信頼できるデータ圧縮が実現できるということだったので、$n$を変化させることで伝送誤差がきちんと小さくなっていくということをシミュレーションにより確認してみます[^6]。

[^6]: ということで、ここで実装する符号化は可逆ではないですし、非可逆としても効率向上の工夫は一切やっていないので、ほぼ実用性のない手法だと思ってください。

今回、前提とした情報源は、$\\{ 0,1 \\}$からなるバイナリ系列です。圧縮法はとても簡単です。伝送したいバイナリ系列全体を$n$個のシンボルからなるブロックに分けます。各ブロックに含まれるバイナリ系列が典型系列かどうかを判定し、典型系列であれば$nR$ビットの符号を典型系列に対して一意になるように割り当てます。そうでなければ諦めて「不明」を表す符号を割り当てます。

### 実装

全体のPythonコードを示します。

```python
import random
import numpy as np
from collections import Counter

def generator(prob, N):

    return [0 if random.random() < prob else 1 for i in range(N)]

def entropy_exp(s):

    prb_0 = list(s).count('0') / len(s)
    prb_1 = 1.0 - prb_0
    if prb_0 == 1.0 or prb_1 == 1.0:
        ent_exp = 0.0
    else:
        ent_exp = - prb_0*np.log2(prb_0)- prb_1*np.log2(prb_1)

    return ent_exp
    
def coder(data, blk, ent, eps, rat):

    # binary data to symbols (= message)
    # ex) [0,1,1,0,1,0,0,1,...] -> ['0110','1001',...] (for blk = 4) 
    message = [''.join(map(str,data[i:i+blk])) for i in range(0,len(data),blk)]

    # histogram
    syms = list(Counter(message).keys())
    frqs = list(Counter(message).values())
    prbs = [frqs[i]/len(message) for i in range(len(frqs))]
    hist = dict(zip(syms, prbs))

    # codebook
    index = 0
    digits = int(rat*blk)  # number of digits for each code
    codebook = {}
    for s,p in hist.items():
        ent_exp = entropy_exp(s)
        if abs(ent_exp - ent) <= eps and index < 2**digits-1:
            codebook[s] = '{:0{digits}b}'.format(index, digits=digits)
            index += 1
    codebook[False] = '{:0{digits}b}'.format(index, digits=digits)

    # coding (message -> codes)
    codes = [codebook[s] if s in codebook else codebook[False]
             for s in message]

    # decodebook
    decodebook = {}
    for k,v in codebook.items():
        decodebook[v] = k
    
    return (codes, decodebook)

def decoder(codes, decodebook, block):

    message = [decodebook[c] for c in codes]
    blocks = [list(s) if s != False else generator(0.5,block) for s in message]
    data_str = sum(blocks, [])

    return [int(d) for d in data_str]

def error_rate(data_in, data_out):

    N = len(data_in)
    count = 0
    for i in range(N):
        if data_in[i] != data_out[i]:
            count += 1
            
    return count / N

if __name__ == "__main__":

    # settings
    BLK = 10          # block size
    NUM = BLK * 1000  # number of binary data sequence
    PRB = 0.8         # probability to generate '0'
    ENT = - PRB * np.log2(PRB) - (1-PRB) * np.log2(1-PRB)  # entropy (= 0.7219..)
    RAT = 0.9         # transmission rate (RAT > ENT)
    EPS = 0.15        # eplsilon for typical sequence
    print("NUM = ", NUM)
    print("BLK = ", BLK)
    print("ENT = ", ENT)
    print("RAT = ", RAT)

    # generate data sequence
    data_in = generator(PRB, NUM)

    # code the data
    (codes, decodebook) = coder(data_in, BLK, ENT, EPS, RAT)

    # decode the data
    data_out = decoder(codes, decodebook, BLK)

    # evaluate error rate
    err_rate = error_rate(data_in, data_out)
    print("error rate = ", err_rate)
```

何をやっているか説明します。main処理部を見てください。

```python
# settings
BLK = 10          # block size
NUM = BLK * 1000  # number of binary data sequence
PRB = 0.8         # probability to generate '0'
ENT = - PRB * np.log2(PRB) - (1-PRB) * np.log2(1-PRB)  # entropy (= 0.7219..)
RAT = 0.9         # transmission rate (RAT > ENT)
EPS = 0.15        # eplsilon for typical sequence
```

ここは各種パラメータを設定する部分です。BLKはブロックのサイズ、すなわち１ブロックに含まれるバイナリ値の数です。NUMは、伝送したい全体のバイナリ値の数です（ブロック1000個分のバイナリ値を伝送する想定にしました）。PRBは"0"が生起する確率です（0.8としました）。ENTは1シンボルあたりのエントロピーです（PRBの値0.8で計算すると、0.7219..になります）。RATは伝送レートです。0.7219..よりも大きくして、かつ、1.0よりも小さくしないとデータ圧縮にならないので、適当な値0.9としました。EPSは典型系列を判定するためのしきい値$\epsilon$の値です。適度に小さい値にしないといけないのですが[^7]、小さくしすぎると、とても大きな$n$を設定しないと誤差が小さくならないので、適当な値として0.15を選びました。

[^7]: $H(X)+\epsilon \leq R$を満たす$\epsilon$にする必要があります。

```python
# generate data sequence
data_in = generator(PRB, NUM)
```

まず、関数generatorで、伝送すべきバイナリ値の系列を作成します。PRBで設定した確率に従って"0"または"1"という整数値をランダムに生成して、長さNUMのリストを出力します。詳細は、関数定義を見てください。

```python
# code the data
(codes, decodebook) = coder(data_in, BLK, ENT, EPS, RAT)
```

関数coderで、data_inを符号化して、符号系列codesと復号化のための辞書decodebookを出力します。関数の中身を見てみます。

```python
# binary data to symbols (= message)
# ex) [0,1,1,0,1,0,0,1,...] -> ['0110','1001',...] (for blk = 4) 
message = [''.join(map(str,data[i:i+blk])) for i in range(0,len(data),blk)]
```

伝送するデータは"0","1"という整数値からなるリストなので、それをブロックごとに文字列にしてmessageというリストにします。ブロックサイズ4の場合の例をコメントに記載しておきました。

```python
# histogram
syms = list(Counter(message).keys())
frqs = list(Counter(message).values())
prbs = [frqs[i]/len(message) for i in range(len(frqs))]
hist = dict(zip(syms, prbs))
```

messageに含まれる要素からヒストグラムを作成し、histという辞書に格納します。

```python
# codebook
index = 0
digits = int(rat*blk)  # number of digits for each code
codebook = {}
for s,p in hist.items():
    ent_exp = entropy_exp(s)
    if abs(ent_exp - ent) <= eps and index < 2**digits-1:
        codebook[s] = '{:0{digits}b}'.format(index, digits=digits)
        index += 1
codebook[False] = '{:0{digits}b}'.format(index, digits=digits)
```

登場するmessageの要素、すなわち、n個のシンボルからなる系列がわかったので、その各々にどんな符号を割り当てるかを決めて、それをcodebookという辞書データに格納します。上記、forループでヒストグラムに含まれているシンボル系列sが典型系列かどうかを判定すべく、関数entropy_expでこの系列に対するエントロピーを計算します。そのエントロピーと真のエントロピーとの差が$\epsilon$以下になっている場合、典型系列とみなし、0.9*nビットの符号indexを0から順に割り当てます。非典型系列の場合、「不明」という意味でFalseとおき、それに対する符号としてindexの最大値を割り当てます。

```python
# coding (message -> codes)
codes = [codebook[s] if s in codebook else codebook[False]
         for s in message]
```

で、得られたコードブックを使って、messageを符号語のリストcodesに変換します。

```python
# decodebook
decodebook = {}
for k,v in codebook.items():
    decodebook[v] = k
```

復号化の際に使えるようにcodebookの逆変換に相当する辞書decodebookを作成します。

```python
return (codes, decodebook)
```

で、codesとdecodebookをリターンします。

main処理部に戻ります。

```python
# decode the data
data_out = decoder(codes, decodebook, BLK)
```

関数decoderで、codesを復号化します。その際、decodebookを使います。関数定義を見ればわかると思いますが、decodebookを参照しながら復号するのですが、復号した結果、Falseだった場合、誤差覚悟でランダムな系列を発生します。

復号化が完了したら、

```python
# evaluate error rate
err_rate = error_rate(data_in, data_out)
print("error rate = ", err_rate)
```

error_rate関数で、エラー率を計算して、表示します。

### 動作確認

では、ブロックサイズ（$n$の値）をいくつか変えてエラー率がどうなるか試してみます。まず、$n=10$の場合です。

```
NUM =  10000
BLK =  10
ENT =  0.7219280948873623
RAT =  0.9
error rate =  0.3499
```

エラー率は、とても大きいです。約35%がエラーです。次に、$n=100$の場合です。

```
NUM =  100000
BLK =  100
ENT =  0.7219280948873623
RAT =  0.9
error rate =  0.03194
```

先程と比べ、１桁減って、3.2%くらいになりました。最後に、$n=200$の場合です。

```
NUM =  200000
BLK =  200
ENT =  0.7219280948873623
RAT =  0.9
error rate =  0.00454
```

さらに１桁減り、0.45%くらいになりました。

ということで、$\epsilon$の値を固定して、$n$の値を大きくしていくと、エラー率がどんどん減っていくことがわかりました。シャノンの定理が主張しているのは、この極限において、エラー率がゼロになるということです。


## おわりに

今回は、古典情報理論のお話だったので、量子計算シミュレータ[qlazy](https://github.com/samn33/qlazy)は使っていません。コードをコピペすれば、お使いのPython環境でも動作すると思いますので、ご興味あれば、遊んでみてください。

次回は、「はじめに」で述べた通り、シャノンの定理の量子版である「雑音のない量子チャネルにおけるシューマッハの符号化定理」を説明します。典型系列に類似した概念として「典型部分空間」というものが出てきますが、今回の議論と同じような流れで、その定理が導かれます。量子情報でも古典と同様に圧縮でき、限界性能も量子エントロピーであるというのは、興味深いです。


以上