Quantum Probabilistic Modelling¶
Quantum computers can be used to calculate joint probability distributions of correlated probabilsistc parameters due to the probabilistic nature. Here we show the idea of how these models can be used for risk modelling.
One example is modelling business risks. An important question for any business is: What would be the impact of external adverse developments on future revenues? For example, these developments could be caused by (macro-)economic, or political events, by competition, or changes in the company’s reputation, or regulation or tax, etc. A major problem for the analysis of this kind of environmental circumstances is that these factors not only (potentially) influence revenues, but are also interlinked in many cases, i.e. they influence each other.
(for further information about business risk see https://arxiv.org/abs/2103.05475)
In [1]:
#
# Import pygrnd and other libraries needed.
# We build on top of the open source framework Qiskit (qiskit.org).
#
import pygrnd
from pygrnd.qc.helper import *
from pygrnd.qc.brm import *
from pygrnd.qc.brm_oracle import *
from pygrnd.qc.QAE import *
from qiskit import execute
from qiskit import Aer
from math import pi
import math
import cmath
import random
from qiskit.visualization import plot_histogram
import matplotlib.pyplot as plt
from scipy.stats import norm
import networkx as nx
from IPython.display import Image
Define and run the model¶
the model is a Graphical Probabilistic Model (GPM), which is a graph with nodes and edges
in its easiest form it is a binary decision tree, equivalent to a directed graph (only one direction) where each node has only one incoming node
extensions allow several edges into one node
risk parameter: nodes=['0','1']
- list of correlations: edges=[('0','1')]
- probabilities that an item gets triggered: probsNodes={'0':0.1,'1':0.1}
- probabilities for transition: probsEdges={('0','1'):0.2}
In [2]:
#
# Small example.
#
nodes=['0','1'] # risk items defition
edges=[('0','1')] # correlations
probsNodes={'0':0.1,'1':0.1} # intrinsic probabilities
probsEdges={('0','1'):0.2} # transition probabilities
print("Number of risk items = ",len(nodes))
print("Number of correlations = ",len(edges))
Number of risk items = 2 Number of correlations = 1
In [3]:
#
# Evaluate the model classically. We obtain a list of probabilities
# and the sum of all probabilities. For two risk items A and B we obtain
# the probabilities for the configurations 00, 01, 10, 11 for (x_B,x_A),
# where x_A and x_B denote the state on/off of A and B respectively. The
# order of the risk items for the list of configurations is inversed to
# be consistent with the Qiskit qubit order (most significant qubit is
# on the bottom of quantum circuits).
#
probs, sumProbs = modelProbabilities(nodes, edges, probsNodes, probsEdges)
print("exact probabilities:",probs)
#
# Convert the list of probabilities to a dictionary to make it
# suitable for comparison with the quantum results below.
#
probsExact={}
for i in range(len(probs)):
cts = format(i, "0"+str(len(nodes))+"b")
probsExact[cts]=probs[i]
exact probabilities: [0.81, 0.072, 0.08999999999999998, 0.027999999999999987]
Monte Carlo simulation of the risk model¶
In [4]:
#
# Compare exact result with three Monte Carlo evaluations of the risk model with 100 rounds each.
#
rounds=100
probsMonteCarlo1=evaluateRiskModelMonteCarlo(nodes, edges, probsNodes, probsEdges, rounds)
probsMonteCarlo2=evaluateRiskModelMonteCarlo(nodes, edges, probsNodes, probsEdges, rounds)
probsMonteCarlo3=evaluateRiskModelMonteCarlo(nodes, edges, probsNodes, probsEdges, rounds)
plot_histogram([probsExact, probsMonteCarlo1, probsMonteCarlo2, probsMonteCarlo3], legend=['exact', 'MC run 1', 'MC run 2', 'MC run 3'], figsize=(16,8))
Out[4]:
In [5]:
#
# Compare exact result with three Monte Carlo evaluations of the risk model with 1000 rounds each.
#
rounds=1000
probsMonteCarlo1=evaluateRiskModelMonteCarlo(nodes, edges, probsNodes, probsEdges, rounds)
probsMonteCarlo2=evaluateRiskModelMonteCarlo(nodes, edges, probsNodes, probsEdges, rounds)
probsMonteCarlo3=evaluateRiskModelMonteCarlo(nodes, edges, probsNodes, probsEdges, rounds)
plot_histogram([probsExact, probsMonteCarlo1, probsMonteCarlo2, probsMonteCarlo3], legend=['exact', 'MC run 1', 'MC run 2', 'MC run 3'], figsize=(16,8))
Out[5]:
Quantum version of the risk model¶
In [6]:
#
# Execute function to build a Qiskit gate with the risk model.
#
rm, mat = brm(nodes, edges, probsNodes, probsEdges, model2gate=False)
#
# Use the statevector simulator to obtain results with high precision.
#
backend = Aer.get_backend('statevector_simulator')
job = execute(rm, backend)
result=job.result()
v=np.asarray(job.result().get_statevector())
#
# Convert statevector to probabilities of measurement results if
# we measure all qubits at the end.
#
probsStatevector={}
for i in range(len(v)):
probsStatevector[format(i, "0"+str(len(nodes))+"b")]=abs(v[i])**2
#
# Use the qasm simulator. This needs a measurement at
# the end to read out the results.
#
rm.measure_all()
display(rm.draw(output='mpl'))
# Use the QASM simulator.
backend = Aer.get_backend('qasm_simulator')
job = execute(rm, backend,shots=10000)
result=job.result()
counts=job.result().get_counts()
plot_histogram([probsExact, probsStatevector, counts], legend=['classical', 'statevector', 'qasm'], figsize=(12,6))
Out[6]:
Visualization of the model¶
In [7]:
#
# Building the graph from the adjacency matrix for visual representation.
#
G = nx.Graph()
G = nx.from_numpy_array(mat)
nx.draw(G, with_labels=True, alpha=0.8, node_size=500)
Quantum Amplitude Estimation¶
In [8]:
#
# Use the oracle to build QAE and evaluate the probabaility for states in 'STATELIST'.
#
name="test"
STATELIST=["11"] # we search for the joint probability where both parameters are triggered
Nshots=1000 # number of shots
QAEqubits=4 # QAE bit precision --> 8 bit strings at the end of the measurement to differentiate 2^8 = 265 bins
rm, mat = brm(nodes, edges, probsNodes, probsEdges, model2gate=True)
ora = brmoracle("ora", rm, len(nodes), 0, STATELIST)
QAE=qae(QAEqubits, len(nodes), len(nodes), rm, ora)
In [9]:
#
# execution with local quantum simulator (qasm simulator)
#
backend = Aer.get_backend('qasm_simulator')
job = execute(QAE, backend,shots=Nshots)
result=job.result()
counts=job.result().get_counts()
plot_histogram(counts,figsize=(16,8))
Out[9]:
In [10]:
#
# Post processing: What is the probability of the states in STATELIST corresponding to the
# two peaks?
#
showQAEoutput(counts, STATELIST, QAEqubits)
Bin with the highest probability: 0001 Number of Bin with the highest probability: 1 The probability of the tail event ['11'] is: 0.03806023374435662
Out[10]:
0.03806023374435662
Bigger example¶
In [11]:
nodes=['0','1','2','3'] # risk items defition
edges=[('0','1'),('0','2'),('1','3'),('2','3')] # risk items correlations
probsNodes={'0':0.1,'1':0.2,'2':0.3,'3':0.4} # intrinsic probs
probsEdges={('0','1'):0.2,('0','2'):0.3,('1','3'):0.4,('2','3'):0.5} # transition probs
In [12]:
#
# Probabilities from exact classical evaluation.
#
probs, sumProbs = modelProbabilities(nodes,edges,probsNodes,probsEdges)
probsExact={}
for i in range(len(probs)):
cts = format(i, "0"+str(len(nodes))+"b")
probsExact[cts]=probs[i]
#
# Probabilities from Monte Carlo simulation with 1000 rounds.
#
rounds=1000
probsMonteCarlo=evaluateRiskModelMonteCarlo(nodes, edges, probsNodes, probsEdges, rounds)
#
# Probabilities from statevector simulation of quantum circuit.
#
rm, mat = brm(nodes, edges, probsNodes, probsEdges, model2gate=False)
rm.draw(output='mpl')
backend = Aer.get_backend('statevector_simulator')
job = execute(rm, backend)
result=job.result()
v=np.asarray(job.result().get_statevector())
probsStatevector={}
for i in range(len(v)):
probsStatevector[format(i, "0"+str(len(nodes))+"b")]=abs(v[i])**2
plot_histogram([probsExact, probsMonteCarlo, probsStatevector], legend=['classical', 'Monte Carlo', 'statevector'], figsize=(16,8))
Out[12]:
