-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathq2.py
140 lines (122 loc) · 3.93 KB
/
q2.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Sat Dec 25 17:10:21 2021
@author: benoitmuller
question 2.
"""
import numpy as np
#import scipy as sp
#import scipy.stats as st
import matplotlib.pyplot as plt
#import time
#import numpy.random as rnd
#from mes_stats import RandomVariable
#from sobol_new import generate_points
from Payoff import Payoff
"""paramètres ou ça marche bien:
np.random.seed(54321) (7,14) (5,6)
"""
# Fix the seed:
np.random.seed(12345)
# File saving options:
save_figures = True # Change the value if needed
if (save_figures==True and
input("Do you really want to save figures"+
" into files?\n(yes/no): ") == "no"):
save_figures = False
alpha=0.01
NN=(2**np.arange(7,12)).astype(int) #(7,14)
mm=(2**np.arange(5,6)).astype(int) #(5,10)
MC1 = np.zeros((2,len(NN)))
MC2 = np.zeros((2,len(NN)))
QMC1 = np.zeros((2,len(NN)))
QMC2 = np.zeros((2,len(NN)))
plt.figure(figsize=(10, 9)) # figsize=(6,4) by default
plt.suptitle("Goal 2: With pre-integration, estimated error" +
" with confidence $1-" + str(alpha) + "$")
for m in mm:
X=Payoff(m)
for j in range(len(NN)):
MC1[:,j], MC2[:,j] = X.PIMC(NN[j],alpha,order=5)
QMC1[:,j], QMC2[:,j] = X.PIQMC(NN[j],alpha,order=5,K=10)
plt.subplot(221)
plt.loglog(NN,MC1[1,:],'.-',label='$m=$'+str(m))
plt.subplot(222)
plt.loglog(NN,QMC1[1,:],'.-',label='$m=$'+str(m))
plt.subplot(223)
plt.loglog(NN,MC2[1,:],'.-',label='$m=$'+str(m))
plt.subplot(224)
plt.loglog(NN,QMC2[1,:],'.-',label='$m=$'+str(m))
# Finest result (biggest m and N):
print("The computed intervals for m =", m,
", N =", NN[-1],"and alpha =",alpha,"are:")
print("V1: PIMC:",MC1[0,-1],"±",MC1[1,-1])
print(" PIQMC:",QMC1[0,-1],"±",QMC1[1,-1])
print("V2: PIMC:",MC2[0,-1],"±",MC2[1,-1])
print(" PIQMC:",QMC2[0,-1],"±",QMC2[1,-1])
# Error plots:
plt.subplot(221)
plt.loglog(NN,NN**(-0.5),"--",label='$1/\sqrt{N}$')
plt.title("Crude Monte Carlo for V1")
plt.xlabel('Sample size N')
plt.ylabel('Estimated error')
plt.legend()
plt.subplot(222)
plt.loglog(NN,NN**(-0.5),"--",label='$1/\sqrt{N}$')
plt.loglog(NN,10/NN,"--",label='$1/N$')
plt.title("Quasi Monte Carlo for V1")
plt.xlabel('Sample size N')
plt.ylabel('Estimated error')
plt.legend()
plt.subplot(223)
plt.loglog(NN,NN**(-0.5),"--",label='$1/\sqrt{N}$')
plt.title("Crude Monte Carlo for V2")
plt.xlabel('Sample size N')
plt.ylabel('Estimated error')
plt.legend()
plt.subplot(224)
plt.loglog(NN,NN**(-0.5),"--",label='$1/\sqrt{N}$')
plt.loglog(NN,10*1/NN,"--",label='$1/N$')
plt.title("Quasi Monte Carlo for V2")
plt.xlabel('Sample size N')
plt.ylabel('Estimated error')
plt.legend()
plt.tight_layout()
if save_figures == True:
plt.savefig('graphics/q2error.pdf')
# The interval plot:
plt.figure(figsize=(10,9))
plt.suptitle("Goal 1: With pre-integration, interval of confidence $1-"
+ str(alpha) + "$, m="+str(mm[-1]))
plt.subplot(221)
plt.xscale('log')
plt.fill_between(NN,+MC1[0,:]+MC1[1,:],MC1[0,:]-MC1[1,:],alpha=0.5)
plt.plot(NN,MC1[0,:],'.-')
plt.title("Crude Monte Carlo for V1")
plt.xlabel('Sample size N')
plt.ylabel('mean and confidence interval')
plt.subplot(222)
plt.xscale('log')
plt.fill_between(NN,QMC1[0,:]+QMC1[1,:],QMC1[0,:]-QMC1[1,:],alpha=0.5)
plt.plot(NN,QMC1[0,:],'.-')
plt.title("Quasi Monte Carlo for V1")
plt.xlabel('Sample size N')
plt.ylabel('mean and confidence interval')
plt.subplot(223)
plt.xscale('log')
plt.fill_between(NN,+MC2[0,:]+MC2[1,:],MC2[0,:]-MC2[1,:],alpha=0.5)
plt.plot(NN,MC2[0,:],'.-')
plt.title("Crude Monte Carlo for V2")
plt.xlabel('Sample size N')
plt.ylabel('mean and confidence interval')
plt.subplot(224)
plt.xscale('log')
plt.fill_between(NN,QMC2[0,:]+QMC2[1,:],QMC2[0,:]-QMC2[1,:],alpha=0.5)
plt.plot(NN,QMC2[0,:],'.-')
plt.title("Quasi Monte Carlo for V2")
plt.xlabel('Sample size N')
plt.ylabel('mean and confidence interval')
plt.tight_layout()
if save_figures == True:
plt.savefig('graphics/q2interval.pdf')