Qiskitで対数正規分布

$$
\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}}
$$
#対数正規分布をqubitで表現
QiskitのAPIはこちら。

確率変数$X$が対数正規分布関数に従っているとき、$X$の確率密度関数は

$$
\mathbb{P}(X=x) = \frac{1}{x\sqrt{2\pi\sigma^2}}\exp\left[-\frac{(\log x - \mu)^2}{\sigma^2}\right]
$$
と表されます。$\mu$は平均、$\sigma$は標準偏差($\sigma^2$が分散)です。

Qiskitではqubitの数を$n$として
$$
\mathcal{P}_X\ket{0}^n = \sum_{i=0}^{2^n-1}\sqrt{\mathbb{P}(x_i)}\ket{i}
$$
を回路に作用することに相当します。qubitでは$\ket{0},,\ket{1}$の二つの状態が観測されるので、$n$個のqubitがあるとき$2^n$個の状態を表現することができます。平方根になっているのは係数の二乗が観測確率であるからです。メソッドの定義は以下の通り。

LogNormalDistribution(num_qubits, mu=None, sigma=None, bounds=None, upto_diag=False, name='P(X)')

• num_qubits: qubitの数
• mu: 平均$\mu$
• sigma: 分散$\sigma^2$(引数の名前はsigmaですが$\sigma^2$の値を入れることに注意)
• bounds: 確率変数の範囲。リストまたはタプル。Noneの場合は(0,1)
• upto_diag: Trueの場合、確率の平方根と対角の乗算を読み込み、回路をより効率的にする
• name: 回路の名前

確率変数の範囲をbounds=[low, high]とするとlowからhighを等間隔に$2^n$個に分割します。muに平均のベクトル、sigmaに共分散行列、boundsに範囲のリストを渡せば多変量対数正規分布を計算することができます。

#IBM Quantumで実験

import matplotlib.pyplot as plt
from qiskit import QuantumCircuit, Aer, IBMQ, transpile, assemble
from qiskit.visualization import plot_histogram
from qiskit_finance.circuit.library import LogNormalDistribution

backend = Aer.get_backend("qasm_simulator")

num_qubits = 3
mu = 3.0
variance = 1.0
low = 0.0
high = 6.0


model = LogNormalDistribution(num_qubits, mu=mu, sigma=variance, bounds=[low, high], name="LogNorm")
x = model.values
y = model.probabilities

# plot probability distribution
plt.bar(x, y, width=(high-low)/2**num_qubits)
plt.xticks(x)
plt.yticks()
plt.grid()
plt.xlabel('random variable', size=15)
plt.ylabel('Probability', size=15)
plt.show()

qubitが3つなのでこの回路では$2^3=8$個の状態を表現することができ、棒が8本見えています(0は観測されない)。また、棒が$\frac{\mbox{high}-\mbox{low}}{8-1}\simeq 0.857$おきになっていることも確認できます。

for i in range(len(x)):
    print("x: {:.3f} y: {:.3f}".format(x[i], y[i]))

x: 0.000 y: 0.000
x: 0.857 y: 0.022
x: 1.714 y: 0.076
x: 2.571 y: 0.126
x: 3.429 y: 0.164
x: 4.286 y: 0.190
x: 5.143 y: 0.206
x: 6.000 y: 0.216

qc = QuantumCircuit(num_qubits, 1)
qc.append(model, [i for i in range(num_qubits)])
#for i in range(num_qubits):
#    qc.measure(i,i)
qc.measure_all()
qc.draw()

t_qc = transpile(qc, backend)
qobj = assemble(t_qc)

counts_100 = backend.run(qobj, shots=100).result().get_counts() 
counts_1k = backend.run(qobj, shots=1000).result().get_counts()
counts_10k = backend.run(qobj, shots=10000).result().get_counts()

plot_histogram(counts_100)

plot_histogram(counts_1k)

plot_histogram(counts_10k)

観測回数(shots)を増やせば、理論値に近づいていくことがわかります。