updating book

This commit is contained in:
Morten Hjorth-Jensen
2021-08-26 10:38:31 +02:00
parent e23be69563
commit c20f947816
66 changed files with 13205 additions and 445 deletions
+36
View File
@@ -0,0 +1,36 @@
Translating doconce text in chapter1.do.txt to ipynb
*** replacing \bm{...} by \boldsymbol{...} (\bm is not supported by MathJax)
*** warning: latex envir \begin{bmatrix} does not work well in Markdown. Stick to \[ ... \], equation, equation*, align, or align* environments in math environments.
*** warning: latex envir \begin{bmatrix} does not work well in Markdown. Stick to \[ ... \], equation, equation*, align, or align* environments in math environments.
Failed to remove ans_at_end environment
Failed to remove sol_at_end environment
output in chapter1.ipynb
Translating doconce text in chapter1.do.txt to ipynb
*** replacing \bm{...} by \boldsymbol{...} (\bm is not supported by MathJax)
*** warning: latex envir \begin{bmatrix} does not work well in Markdown. Stick to \[ ... \], equation, equation*, align, or align* environments in math environments.
*** warning: latex envir \begin{bmatrix} does not work well in Markdown. Stick to \[ ... \], equation, equation*, align, or align* environments in math environments.
Failed to remove ans_at_end environment
Failed to remove sol_at_end environment
output in chapter1.ipynb
Translating doconce text in chapter1.do.txt to ipynb
*** replacing \bm{...} by \boldsymbol{...} (\bm is not supported by MathJax)
*** warning: latex envir \begin{bmatrix} does not work well in Markdown. Stick to \[ ... \], equation, equation*, align, or align* environments in math environments.
*** warning: latex envir \begin{bmatrix} does not work well in Markdown. Stick to \[ ... \], equation, equation*, align, or align* environments in math environments.
Failed to remove ans_at_end environment
Failed to remove sol_at_end environment
output in chapter1.ipynb
Translating doconce text in chapter1.do.txt to ipynb
*** replacing \bm{...} by \boldsymbol{...} (\bm is not supported by MathJax)
*** warning: latex envir \begin{bmatrix} does not work well in Markdown. Stick to \[ ... \], equation, equation*, align, or align* environments in math environments.
*** warning: latex envir \begin{bmatrix} does not work well in Markdown. Stick to \[ ... \], equation, equation*, align, or align* environments in math environments.
Failed to remove ans_at_end environment
Failed to remove sol_at_end environment
output in chapter1.ipynb
+5 -52
View File
@@ -1263,22 +1263,22 @@ matrices as upper case boldfaced letters.
!bt
\[
\frac{\partial (\bm{b}^T\bm{a})}{\partial \bm{a}} = \bm{b},
\frac{\partial (\bm{b}^T\bm{a})}{\partial\bm{a}}=\bm{b},
\]
!et
!bt
\[
\frac{\partial (\bm{a}^T\bm{A}\bm{a})}{\partial \bm{a}} = (\bm{A}+\bm{A}^T)\bm{a},
\frac{\partial (\bm{a}^T\bm{A}\bm{a})}{\partial\bm{a}}=(\bm{A}+\bm{A}^T)\bm{a},
\]
!et
!bt
\[
\frac{\partial tr(\bm{B}\bm{A})}{\partial \bm{A}} = \bm{B}^T,
\frac{\partial tr(\bm{B}\bm{A})}{\partial\bm{A}}=\bm{B}^T,
\]
!et
!bt
\[
\frac{\partial \log{\vert\bm{A}\vert}}{\partial \bm{A}} = (\bm{A}^{-1})^T.
\frac{\partial\log{\vert\bm{A}\vert}}{\partial \bm{A}}=(\bm{A}^{-1})^T.
\]
!et
@@ -2137,51 +2137,6 @@ where we have defined the mean value of $\hat{y}$ as
You can use the functionality included in scikit-learn. If you feel for it, you can use your own program and define functions which compute the above two functions.
Discuss the meaning of these results. Try also to vary the coefficient in front of the added stochastic noise term and discuss the quality of the fits.
!bsol
The code here is an example of where we define our own design matrix and fit parameters $\beta$.
!bc pycod
import os
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split
def save_fig(fig_id):
plt.savefig(image_path(fig_id) + ".png", format='png')
def R2(y_data, y_model):
return 1 - np.sum((y_data - y_model) ** 2) / np.sum((y_data - np.mean(y_data)) ** 2)
def MSE(y_data,y_model):
n = np.size(y_model)
return np.sum((y_data-y_model)**2)/n
x = np.random.rand(100)
y = 2.0+5*x*x+0.1*np.random.randn(100)
# The design matrix now as function of a given polynomial
X = np.zeros((len(x),3))
X[:,0] = 1.0
X[:,1] = x
X[:,2] = x**2
# We split the data in test and training data
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)
# matrix inversion to find beta
beta = np.linalg.inv(X_train.T @ X_train) @ X_train.T @ y_train
print(beta)
# and then make the prediction
ytilde = X_train @ beta
print("Training R2")
print(R2(y_train,ytilde))
print("Training MSE")
print(MSE(y_train,ytilde))
ypredict = X_test @ beta
print("Test R2")
print(R2(y_test,ypredict))
print("Test MSE")
print(MSE(y_test,ypredict))
!ec
!esol
@@ -2266,8 +2221,6 @@ Perform an ordinary least squares and compute the means squared error and the $R
!bsubex
Add now a model which allows you to make polynomials up to degree $15$. Perform a standard OLS fitting of the training data and compute the MSE and $R2$ for the training and test data and plot both test and training data MSE and $R2$ as functions of the polynomial degree. Compare what you see with Figure 2.11 of Hastie et al. Comment your results. For which polynomial degree do you find an optimal MSE (smallest value)?
!bsol
Here you simply need to change the degree of the polynomial in the above code to $n=15$.
!esol
!esubex
+90
View File
@@ -0,0 +1,90 @@
3.3773726001100143E-005, 3.1715032621225665E-002
2.7018980800880114E-004, 0.25379346864959029
9.1189060202970370E-004, 0.85691856357595775
2.1615184640704091E-003, 2.0322700794292126
4.2217157501375172E-003, 3.9716946611295496
7.2951248162376296E-003, 6.8678138410008760
1.1584388018377344E-002, 10.913973030767592
1.7292147712563273E-002, 16.304871533080519
2.4621046254802003E-002, 23.235226466525493
3.3773726001100138E-002, 31.902701432416372
4.4952829307464297E-002, 42.504284679798289
5.8360998529901037E-002, 55.243642496340065
7.4200876024417023E-002, 70.325036648868448
9.2675104147018753E-002, 87.967299710326188
0.11398632525371297, 108.39264632432710
0.13833718170050618 , 131.84282952216495
0.16593031584340495 , 158.57124952965228
0.19696837003841602 , 188.82820513039681
0.203199995458126059, 195.042764094891368
0.211199995279312130, 202.937172130123315
0.219199995100498229, 210.852580165355278
0.227199994921684245, 218.789988200587175
0.235199994742870316, 226.748396235819115
0.243199994564056388, 234.740804271051104
0.251199994385242487, 242.757212306283037
0.259199994206428530, 250.797620341515000
0.267199994027614574, 258.874028376746878
0.275199993848800673, 266.977436411978829
0.283199993669986716, 275.092844447210780
0.291199993491172815, 283.248252482442751
0.299199993312358858, 291.445660517674696
0.307199993133544902, 299.655068552906641
0.315199992954731001, 307.909476588138546
0.323199992775917044, 316.192884623370503
0.339199992418289187, 332.846700693834407
0.355199992060661329, 349.620516764298316
0.371199991703033416, 366.537332834762140
0.387199991345405559, 383.604148905225998
0.403199990987777701, 400.801964975689941
0.419199990630149844, 418.158781046153820
0.435199990272521986, 435.624597116617792
0.451199989914894073, 453.288413187081574
0.467199989557266215, 471.062229257545425
0.483199989199638358, 489.048045328009380
0.499199988842010500, 507.146861398473277
0.515199988484382643, 525.392677468937222
0.531199988126754730, 543.829493539401028
0.531999988108873390, 544.796634342924222
0.550079987704753859, 565.860416502548446
0.567999987304210641, 586.865970501467928
0.585599986910820047, 607.703068178978242
0.603199986517429343, 628.645165856488575
0.620639986127614951, 649.578035373294142
0.637919985741376872, 670.496676729395176
0.655199985355138792, 691.549318085496111
0.672479984968900713, 712.766959441597237
0.689599984586238834, 733.981372636993456
0.706879984200000755, 755.521013993094471
0.724159983813762675, 777.228655349195492
0.741439983427524596, 799.117296705296553
0.758719983041286516, 821.191938061397536
0.775999982655048326, 843.460579417498366
0.793439982265233934, 866.080448934304059
0.810879981875419542, 888.907318451109631
0.828479981482028949, 912.100416128620054
0.846239981085061932, 935.664741966834868
0.864159980684518825, 959.607295965754474
0.882079980283975607, 983.787849964674024
0.900159979879856187, 1008.36063212429838
0.918399979472160344, 1033.33564244462718
0.936639979064464612, 1058.56865276495591
0.955199978649616255, 1084.36811940669395
0.973919978231191585, 1110.59381420913678
0.992799977809190715, 1137.25073717228429
1.01183997738361353, 1164.35288829613614
1.06047997629642499, 1234.41724915034661
1.11039997518062594, 1307.72443529019392
1.16191997402906422, 1384.74890303708753
1.21535997283458719, 1466.00110871243692
1.27071997159719463, 1551.73305231624181
1.32847997030615828, 1642.68741833061654
1.38895996895432461, 1739.56266307697001
1.45263996753096580, 1843.32947103741640
1.52031996601820008, 1955.19398301547858
1.59295996439456933, 2077.14256797538474
1.67167996263504026, 2211.44382304206692
1.75855996069312082, 2361.74471430468566
1.85647995850443825, 2533.64534865592441
1.96959995597600934, 2734.93465827410455
2.10543995293974895, 2980.03336671234320
1 3.3773726001100143E-005 3.1715032621225665E-002
2 2.7018980800880114E-004 0.25379346864959029
3 9.1189060202970370E-004 0.85691856357595775
4 2.1615184640704091E-003 2.0322700794292126
5 4.2217157501375172E-003 3.9716946611295496
6 7.2951248162376296E-003 6.8678138410008760
7 1.1584388018377344E-002 10.913973030767592
8 1.7292147712563273E-002 16.304871533080519
9 2.4621046254802003E-002 23.235226466525493
10 3.3773726001100138E-002 31.902701432416372
11 4.4952829307464297E-002 42.504284679798289
12 5.8360998529901037E-002 55.243642496340065
13 7.4200876024417023E-002 70.325036648868448
14 9.2675104147018753E-002 87.967299710326188
15 0.11398632525371297 108.39264632432710
16 0.13833718170050618 131.84282952216495
17 0.16593031584340495 158.57124952965228
18 0.19696837003841602 188.82820513039681
19 0.203199995458126059 195.042764094891368
20 0.211199995279312130 202.937172130123315
21 0.219199995100498229 210.852580165355278
22 0.227199994921684245 218.789988200587175
23 0.235199994742870316 226.748396235819115
24 0.243199994564056388 234.740804271051104
25 0.251199994385242487 242.757212306283037
26 0.259199994206428530 250.797620341515000
27 0.267199994027614574 258.874028376746878
28 0.275199993848800673 266.977436411978829
29 0.283199993669986716 275.092844447210780
30 0.291199993491172815 283.248252482442751
31 0.299199993312358858 291.445660517674696
32 0.307199993133544902 299.655068552906641
33 0.315199992954731001 307.909476588138546
34 0.323199992775917044 316.192884623370503
35 0.339199992418289187 332.846700693834407
36 0.355199992060661329 349.620516764298316
37 0.371199991703033416 366.537332834762140
38 0.387199991345405559 383.604148905225998
39 0.403199990987777701 400.801964975689941
40 0.419199990630149844 418.158781046153820
41 0.435199990272521986 435.624597116617792
42 0.451199989914894073 453.288413187081574
43 0.467199989557266215 471.062229257545425
44 0.483199989199638358 489.048045328009380
45 0.499199988842010500 507.146861398473277
46 0.515199988484382643 525.392677468937222
47 0.531199988126754730 543.829493539401028
48 0.531999988108873390 544.796634342924222
49 0.550079987704753859 565.860416502548446
50 0.567999987304210641 586.865970501467928
51 0.585599986910820047 607.703068178978242
52 0.603199986517429343 628.645165856488575
53 0.620639986127614951 649.578035373294142
54 0.637919985741376872 670.496676729395176
55 0.655199985355138792 691.549318085496111
56 0.672479984968900713 712.766959441597237
57 0.689599984586238834 733.981372636993456
58 0.706879984200000755 755.521013993094471
59 0.724159983813762675 777.228655349195492
60 0.741439983427524596 799.117296705296553
61 0.758719983041286516 821.191938061397536
62 0.775999982655048326 843.460579417498366
63 0.793439982265233934 866.080448934304059
64 0.810879981875419542 888.907318451109631
65 0.828479981482028949 912.100416128620054
66 0.846239981085061932 935.664741966834868
67 0.864159980684518825 959.607295965754474
68 0.882079980283975607 983.787849964674024
69 0.900159979879856187 1008.36063212429838
70 0.918399979472160344 1033.33564244462718
71 0.936639979064464612 1058.56865276495591
72 0.955199978649616255 1084.36811940669395
73 0.973919978231191585 1110.59381420913678
74 0.992799977809190715 1137.25073717228429
75 1.01183997738361353 1164.35288829613614
76 1.06047997629642499 1234.41724915034661
77 1.11039997518062594 1307.72443529019392
78 1.16191997402906422 1384.74890303708753
79 1.21535997283458719 1466.00110871243692
80 1.27071997159719463 1551.73305231624181
81 1.32847997030615828 1642.68741833061654
82 1.38895996895432461 1739.56266307697001
83 1.45263996753096580 1843.32947103741640
84 1.52031996601820008 1955.19398301547858
85 1.59295996439456933 2077.14256797538474
86 1.67167996263504026 2211.44382304206692
87 1.75855996069312082 2361.74471430468566
88 1.85647995850443825 2533.64534865592441
89 1.96959995597600934 2734.93465827410455
90 2.10543995293974895 2980.03336671234320
File diff suppressed because it is too large Load Diff
File diff suppressed because it is too large Load Diff
@@ -0,0 +1,67 @@
// This function computes the autocorrelation function for
// the standard c++ random number generator
#include <fstream>
#include <iomanip>
#include <iostream>
#include <cmath>
using namespace std;
// output file as global variable
ofstream ofile;
// Main function begins here
int main(int argc, char* argv[])
{
int n;
char *outfilename;
cin >> n;
double MCint = 0.; double MCintsqr2=0.;
double invers_period = 1./RAND_MAX; // initialise the random number generator
srand(time(NULL)); // This produces the so-called seed in MC jargon
// Compute the variance and the mean value of the uniform distribution
// Compute also the specific values x for each cycle in order to be able to
// the covariance and the correlation function
// Read in output file, abort if there are too few command-line arguments
if( argc <= 2 ){
cout << "Bad Usage: " << argv[0] <<
" read also output file and number of cycles on same line" << endl;
exit(1);
}
else{
outfilename=argv[1];
}
ofile.open(outfilename);
// Get the number of Monte-Carlo samples
n = atoi(argv[2]);
double *X;
X = new double[n];
for (int i = 0; i < n; i++){
double x = double(rand())*invers_period;
X[i] = x;
MCint += x;
MCintsqr2 += x*x;
}
double Mean = MCint/((double) n );
MCintsqr2 = MCintsqr2/((double) n );
double STDev = sqrt(MCintsqr2-Mean*Mean);
double Variance = MCintsqr2-Mean*Mean;
// Write mean value and standard deviation
cout << " Standard deviation= " << STDev << " Integral = " << Mean << endl;
// Now we compute the autocorrelation function, setting the distance d between two
// to a most 1/4 of the total number of cycles
double *autocor; autocor = new double[n];
for (int j = 0; j < n; j++){
double sum = 0.0;
for (int k = 0; k < (n-j); k++){
sum += (X[k]-Mean)*(X[k+j]-Mean);
}
autocor[j] = sum/Variance/((double) n );
ofile << setiosflags(ios::showpoint | ios::uppercase);
ofile << setw(15) << setprecision(8) << j;
ofile << setw(15) << setprecision(8) << autocor[j] << endl;
}
ofile.close(); // close output file
return 0;
} // end of main program
@@ -0,0 +1,76 @@
// This function computes the autocorrelation function for
// the standard c++ random number generator
#include <fstream>
#include <iomanip>
#include <iostream>
#include <cmath>
#include <random>
using namespace std;
// output file as global variable
ofstream ofile;
// Main function begins here
int main(int argc, char* argv[])
{
int n;
char *outfilename;
cin >> n;
double MCint = 0.; double MCintsqr2=0.;
// Initialize the seed and call the Mersienne algo
std::random_device rd;
std::mt19937_64 gen(rd());
// Set up the uniform distribution for x \in [[0, 1]
std::uniform_real_distribution<double> RandomNumberGenerator(0.0,1.0);
// Compute the variance and the mean value of the uniform distribution
// Compute also the specific values x for each cycle in order to be able to
// the covariance and the correlation function
// Read in output file, abort if there are too few command-line arguments
if( argc <= 2 ){
cout << "Bad Usage: " << argv[0] <<
" read also output file and number of cycles on same line" << endl;
exit(1);
}
else{
outfilename=argv[1];
}
ofile.open(outfilename);
// Get the number of Monte-Carlo samples
n = atoi(argv[2]);
double *X;
X = new double[n];
for (int i = 0; i < n; i++){
double x = RandomNumberGenerator(gen);
X[i] = x;
MCint += x;
MCintsqr2 += x*x;
}
double Mean = MCint/((double) n );
MCintsqr2 = MCintsqr2/((double) n );
double STDev = sqrt(MCintsqr2-Mean*Mean);
double Variance = MCintsqr2-Mean*Mean;
// Write mean value and standard deviation
cout << " Standard deviation= " << STDev << " Integral = " << Mean << endl;
// Now we compute the autocorrelation function, setting the distance d
double *autocor; autocor = new double[n];
for (int j = 0; j < n; j++){
double sum = 0.0;
for (int k = 0; k < (n-j); k++){
sum += (X[k]-Mean)*(X[k+j]-Mean);
}
autocor[j] = sum/Variance/((double) n );
ofile << setiosflags(ios::showpoint | ios::uppercase);
ofile << setw(15) << setprecision(8) << j;
ofile << setw(15) << setprecision(8) << autocor[j] << endl;
}
ofile.close(); // close output file
return 0;
} // end of main program
File diff suppressed because it is too large Load Diff
+16
View File
@@ -0,0 +1,16 @@
import numpy as np
from matplotlib import pyplot as plt
# Load in data file
data = np.loadtxt("autocor.dat")
data1 = np.loadtxt("automersenne.dat")
# Make arrays containing x-axis and binding energies as function of A
x = data[:,0]
corr = data[:,1]
corr2 = data1[:,1]
plt.plot(x, corr ,'ro', x, corr2, 'b')
plt.axis([0,1000,-0.2, 1.1])
plt.xlabel(r'$d$')
plt.ylabel(r'$C_d$')
plt.title(r'autocorrelation function for RNG')
plt.savefig('autocorr.pdf')
plt.show()
+90
View File
@@ -0,0 +1,90 @@
3.3773726001100143E-005, 3.1715032621225665E-002
2.7018980800880114E-004, 0.25379346864959029
9.1189060202970370E-004, 0.85691856357595775
2.1615184640704091E-003, 2.0322700794292126
4.2217157501375172E-003, 3.9716946611295496
7.2951248162376296E-003, 6.8678138410008760
1.1584388018377344E-002, 10.913973030767592
1.7292147712563273E-002, 16.304871533080519
2.4621046254802003E-002, 23.235226466525493
3.3773726001100138E-002, 31.902701432416372
4.4952829307464297E-002, 42.504284679798289
5.8360998529901037E-002, 55.243642496340065
7.4200876024417023E-002, 70.325036648868448
9.2675104147018753E-002, 87.967299710326188
0.11398632525371297, 108.39264632432710
0.13833718170050618 , 131.84282952216495
0.16593031584340495 , 158.57124952965228
0.19696837003841602 , 188.82820513039681
0.203199995458126059, 195.042764094891368
0.211199995279312130, 202.937172130123315
0.219199995100498229, 210.852580165355278
0.227199994921684245, 218.789988200587175
0.235199994742870316, 226.748396235819115
0.243199994564056388, 234.740804271051104
0.251199994385242487, 242.757212306283037
0.259199994206428530, 250.797620341515000
0.267199994027614574, 258.874028376746878
0.275199993848800673, 266.977436411978829
0.283199993669986716, 275.092844447210780
0.291199993491172815, 283.248252482442751
0.299199993312358858, 291.445660517674696
0.307199993133544902, 299.655068552906641
0.315199992954731001, 307.909476588138546
0.323199992775917044, 316.192884623370503
0.339199992418289187, 332.846700693834407
0.355199992060661329, 349.620516764298316
0.371199991703033416, 366.537332834762140
0.387199991345405559, 383.604148905225998
0.403199990987777701, 400.801964975689941
0.419199990630149844, 418.158781046153820
0.435199990272521986, 435.624597116617792
0.451199989914894073, 453.288413187081574
0.467199989557266215, 471.062229257545425
0.483199989199638358, 489.048045328009380
0.499199988842010500, 507.146861398473277
0.515199988484382643, 525.392677468937222
0.531199988126754730, 543.829493539401028
0.531999988108873390, 544.796634342924222
0.550079987704753859, 565.860416502548446
0.567999987304210641, 586.865970501467928
0.585599986910820047, 607.703068178978242
0.603199986517429343, 628.645165856488575
0.620639986127614951, 649.578035373294142
0.637919985741376872, 670.496676729395176
0.655199985355138792, 691.549318085496111
0.672479984968900713, 712.766959441597237
0.689599984586238834, 733.981372636993456
0.706879984200000755, 755.521013993094471
0.724159983813762675, 777.228655349195492
0.741439983427524596, 799.117296705296553
0.758719983041286516, 821.191938061397536
0.775999982655048326, 843.460579417498366
0.793439982265233934, 866.080448934304059
0.810879981875419542, 888.907318451109631
0.828479981482028949, 912.100416128620054
0.846239981085061932, 935.664741966834868
0.864159980684518825, 959.607295965754474
0.882079980283975607, 983.787849964674024
0.900159979879856187, 1008.36063212429838
0.918399979472160344, 1033.33564244462718
0.936639979064464612, 1058.56865276495591
0.955199978649616255, 1084.36811940669395
0.973919978231191585, 1110.59381420913678
0.992799977809190715, 1137.25073717228429
1.01183997738361353, 1164.35288829613614
1.06047997629642499, 1234.41724915034661
1.11039997518062594, 1307.72443529019392
1.16191997402906422, 1384.74890303708753
1.21535997283458719, 1466.00110871243692
1.27071997159719463, 1551.73305231624181
1.32847997030615828, 1642.68741833061654
1.38895996895432461, 1739.56266307697001
1.45263996753096580, 1843.32947103741640
1.52031996601820008, 1955.19398301547858
1.59295996439456933, 2077.14256797538474
1.67167996263504026, 2211.44382304206692
1.75855996069312082, 2361.74471430468566
1.85647995850443825, 2533.64534865592441
1.96959995597600934, 2734.93465827410455
2.10543995293974895, 2980.03336671234320
1 3.3773726001100143E-005 3.1715032621225665E-002
2 2.7018980800880114E-004 0.25379346864959029
3 9.1189060202970370E-004 0.85691856357595775
4 2.1615184640704091E-003 2.0322700794292126
5 4.2217157501375172E-003 3.9716946611295496
6 7.2951248162376296E-003 6.8678138410008760
7 1.1584388018377344E-002 10.913973030767592
8 1.7292147712563273E-002 16.304871533080519
9 2.4621046254802003E-002 23.235226466525493
10 3.3773726001100138E-002 31.902701432416372
11 4.4952829307464297E-002 42.504284679798289
12 5.8360998529901037E-002 55.243642496340065
13 7.4200876024417023E-002 70.325036648868448
14 9.2675104147018753E-002 87.967299710326188
15 0.11398632525371297 108.39264632432710
16 0.13833718170050618 131.84282952216495
17 0.16593031584340495 158.57124952965228
18 0.19696837003841602 188.82820513039681
19 0.203199995458126059 195.042764094891368
20 0.211199995279312130 202.937172130123315
21 0.219199995100498229 210.852580165355278
22 0.227199994921684245 218.789988200587175
23 0.235199994742870316 226.748396235819115
24 0.243199994564056388 234.740804271051104
25 0.251199994385242487 242.757212306283037
26 0.259199994206428530 250.797620341515000
27 0.267199994027614574 258.874028376746878
28 0.275199993848800673 266.977436411978829
29 0.283199993669986716 275.092844447210780
30 0.291199993491172815 283.248252482442751
31 0.299199993312358858 291.445660517674696
32 0.307199993133544902 299.655068552906641
33 0.315199992954731001 307.909476588138546
34 0.323199992775917044 316.192884623370503
35 0.339199992418289187 332.846700693834407
36 0.355199992060661329 349.620516764298316
37 0.371199991703033416 366.537332834762140
38 0.387199991345405559 383.604148905225998
39 0.403199990987777701 400.801964975689941
40 0.419199990630149844 418.158781046153820
41 0.435199990272521986 435.624597116617792
42 0.451199989914894073 453.288413187081574
43 0.467199989557266215 471.062229257545425
44 0.483199989199638358 489.048045328009380
45 0.499199988842010500 507.146861398473277
46 0.515199988484382643 525.392677468937222
47 0.531199988126754730 543.829493539401028
48 0.531999988108873390 544.796634342924222
49 0.550079987704753859 565.860416502548446
50 0.567999987304210641 586.865970501467928
51 0.585599986910820047 607.703068178978242
52 0.603199986517429343 628.645165856488575
53 0.620639986127614951 649.578035373294142
54 0.637919985741376872 670.496676729395176
55 0.655199985355138792 691.549318085496111
56 0.672479984968900713 712.766959441597237
57 0.689599984586238834 733.981372636993456
58 0.706879984200000755 755.521013993094471
59 0.724159983813762675 777.228655349195492
60 0.741439983427524596 799.117296705296553
61 0.758719983041286516 821.191938061397536
62 0.775999982655048326 843.460579417498366
63 0.793439982265233934 866.080448934304059
64 0.810879981875419542 888.907318451109631
65 0.828479981482028949 912.100416128620054
66 0.846239981085061932 935.664741966834868
67 0.864159980684518825 959.607295965754474
68 0.882079980283975607 983.787849964674024
69 0.900159979879856187 1008.36063212429838
70 0.918399979472160344 1033.33564244462718
71 0.936639979064464612 1058.56865276495591
72 0.955199978649616255 1084.36811940669395
73 0.973919978231191585 1110.59381420913678
74 0.992799977809190715 1137.25073717228429
75 1.01183997738361353 1164.35288829613614
76 1.06047997629642499 1234.41724915034661
77 1.11039997518062594 1307.72443529019392
78 1.16191997402906422 1384.74890303708753
79 1.21535997283458719 1466.00110871243692
80 1.27071997159719463 1551.73305231624181
81 1.32847997030615828 1642.68741833061654
82 1.38895996895432461 1739.56266307697001
83 1.45263996753096580 1843.32947103741640
84 1.52031996601820008 1955.19398301547858
85 1.59295996439456933 2077.14256797538474
86 1.67167996263504026 2211.44382304206692
87 1.75855996069312082 2361.74471430468566
88 1.85647995850443825 2533.64534865592441
89 1.96959995597600934 2734.93465827410455
90 2.10543995293974895 2980.03336671234320
File diff suppressed because it is too large Load Diff
File diff suppressed because it is too large Load Diff
@@ -0,0 +1,67 @@
// This function computes the autocorrelation function for
// the standard c++ random number generator
#include <fstream>
#include <iomanip>
#include <iostream>
#include <cmath>
using namespace std;
// output file as global variable
ofstream ofile;
// Main function begins here
int main(int argc, char* argv[])
{
int n;
char *outfilename;
cin >> n;
double MCint = 0.; double MCintsqr2=0.;
double invers_period = 1./RAND_MAX; // initialise the random number generator
srand(time(NULL)); // This produces the so-called seed in MC jargon
// Compute the variance and the mean value of the uniform distribution
// Compute also the specific values x for each cycle in order to be able to
// the covariance and the correlation function
// Read in output file, abort if there are too few command-line arguments
if( argc <= 2 ){
cout << "Bad Usage: " << argv[0] <<
" read also output file and number of cycles on same line" << endl;
exit(1);
}
else{
outfilename=argv[1];
}
ofile.open(outfilename);
// Get the number of Monte-Carlo samples
n = atoi(argv[2]);
double *X;
X = new double[n];
for (int i = 0; i < n; i++){
double x = double(rand())*invers_period;
X[i] = x;
MCint += x;
MCintsqr2 += x*x;
}
double Mean = MCint/((double) n );
MCintsqr2 = MCintsqr2/((double) n );
double STDev = sqrt(MCintsqr2-Mean*Mean);
double Variance = MCintsqr2-Mean*Mean;
// Write mean value and standard deviation
cout << " Standard deviation= " << STDev << " Integral = " << Mean << endl;
// Now we compute the autocorrelation function, setting the distance d between two
// to a most 1/4 of the total number of cycles
double *autocor; autocor = new double[n];
for (int j = 0; j < n; j++){
double sum = 0.0;
for (int k = 0; k < (n-j); k++){
sum += (X[k]-Mean)*(X[k+j]-Mean);
}
autocor[j] = sum/Variance/((double) n );
ofile << setiosflags(ios::showpoint | ios::uppercase);
ofile << setw(15) << setprecision(8) << j;
ofile << setw(15) << setprecision(8) << autocor[j] << endl;
}
ofile.close(); // close output file
return 0;
} // end of main program
@@ -0,0 +1,76 @@
// This function computes the autocorrelation function for
// the standard c++ random number generator
#include <fstream>
#include <iomanip>
#include <iostream>
#include <cmath>
#include <random>
using namespace std;
// output file as global variable
ofstream ofile;
// Main function begins here
int main(int argc, char* argv[])
{
int n;
char *outfilename;
cin >> n;
double MCint = 0.; double MCintsqr2=0.;
// Initialize the seed and call the Mersienne algo
std::random_device rd;
std::mt19937_64 gen(rd());
// Set up the uniform distribution for x \in [[0, 1]
std::uniform_real_distribution<double> RandomNumberGenerator(0.0,1.0);
// Compute the variance and the mean value of the uniform distribution
// Compute also the specific values x for each cycle in order to be able to
// the covariance and the correlation function
// Read in output file, abort if there are too few command-line arguments
if( argc <= 2 ){
cout << "Bad Usage: " << argv[0] <<
" read also output file and number of cycles on same line" << endl;
exit(1);
}
else{
outfilename=argv[1];
}
ofile.open(outfilename);
// Get the number of Monte-Carlo samples
n = atoi(argv[2]);
double *X;
X = new double[n];
for (int i = 0; i < n; i++){
double x = RandomNumberGenerator(gen);
X[i] = x;
MCint += x;
MCintsqr2 += x*x;
}
double Mean = MCint/((double) n );
MCintsqr2 = MCintsqr2/((double) n );
double STDev = sqrt(MCintsqr2-Mean*Mean);
double Variance = MCintsqr2-Mean*Mean;
// Write mean value and standard deviation
cout << " Standard deviation= " << STDev << " Integral = " << Mean << endl;
// Now we compute the autocorrelation function, setting the distance d
double *autocor; autocor = new double[n];
for (int j = 0; j < n; j++){
double sum = 0.0;
for (int k = 0; k < (n-j); k++){
sum += (X[k]-Mean)*(X[k+j]-Mean);
}
autocor[j] = sum/Variance/((double) n );
ofile << setiosflags(ios::showpoint | ios::uppercase);
ofile << setw(15) << setprecision(8) << j;
ofile << setw(15) << setprecision(8) << autocor[j] << endl;
}
ofile.close(); // close output file
return 0;
} // end of main program
File diff suppressed because it is too large Load Diff
+16
View File
@@ -0,0 +1,16 @@
import numpy as np
from matplotlib import pyplot as plt
# Load in data file
data = np.loadtxt("autocor.dat")
data1 = np.loadtxt("automersenne.dat")
# Make arrays containing x-axis and binding energies as function of A
x = data[:,0]
corr = data[:,1]
corr2 = data1[:,1]
plt.plot(x, corr ,'ro', x, corr2, 'b')
plt.axis([0,1000,-0.2, 1.1])
plt.xlabel(r'$d$')
plt.ylabel(r'$C_d$')
plt.title(r'autocorrelation function for RNG')
plt.savefig('autocorr.pdf')
plt.show()
Binary file not shown.

After

Width:  |  Height:  |  Size: 22 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 16 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 16 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 21 KiB

Binary file not shown.
Binary file not shown.
Binary file not shown.

Before

Width:  |  Height:  |  Size: 6.1 KiB

After

Width:  |  Height:  |  Size: 6.2 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 16 KiB

Binary file not shown.

Before

Width:  |  Height:  |  Size: 9.3 KiB

After

Width:  |  Height:  |  Size: 9.4 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 21 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 25 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 95 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 33 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 10 KiB

Binary file not shown.

Before

Width:  |  Height:  |  Size: 14 KiB

After

Width:  |  Height:  |  Size: 13 KiB

Binary file not shown.

Before

Width:  |  Height:  |  Size: 9.8 KiB

After

Width:  |  Height:  |  Size: 9.9 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 16 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 19 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 26 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 25 KiB

@@ -1963,7 +1963,7 @@
"metadata": {},
"source": [
"$$\n",
"\\frac{\\partial \\log{\\vert\\boldsymbol{A}\\vert}}{\\partial \\boldsymbol{A}} = (\\boldsymbol{A}^{-1})^T.\n",
"\\frac{\\partial\\log{\\vert\\boldsymbol{A}\\vert}}{\\partial \\boldsymbol{A}}=(\\boldsymbol{A}^{-1})^T.\n",
"$$"
]
},
@@ -3359,66 +3359,6 @@
"You can use the functionality included in scikit-learn. If you feel for it, you can use your own program and define functions which compute the above two functions. \n",
"Discuss the meaning of these results. Try also to vary the coefficient in front of the added stochastic noise term and discuss the quality of the fits.\n",
"\n",
"!bsol\n",
"The code here is an example of where we define our own design matrix and fit parameters $\\beta$."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false,
"editable": true
},
"outputs": [],
"source": [
"import os\n",
"import numpy as np\n",
"import pandas as pd\n",
"import matplotlib.pyplot as plt\n",
"from sklearn.model_selection import train_test_split\n",
"\n",
"def save_fig(fig_id):\n",
" plt.savefig(image_path(fig_id) + \".png\", format='png')\n",
"\n",
"def R2(y_data, y_model):\n",
" return 1 - np.sum((y_data - y_model) ** 2) / np.sum((y_data - np.mean(y_data)) ** 2)\n",
"def MSE(y_data,y_model):\n",
" n = np.size(y_model)\n",
" return np.sum((y_data-y_model)**2)/n\n",
"\n",
"x = np.random.rand(100)\n",
"y = 2.0+5*x*x+0.1*np.random.randn(100)\n",
"\n",
"\n",
"# The design matrix now as function of a given polynomial\n",
"X = np.zeros((len(x),3))\n",
"X[:,0] = 1.0\n",
"X[:,1] = x\n",
"X[:,2] = x**2\n",
"# We split the data in test and training data\n",
"X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)\n",
"# matrix inversion to find beta\n",
"beta = np.linalg.inv(X_train.T @ X_train) @ X_train.T @ y_train\n",
"print(beta)\n",
"# and then make the prediction\n",
"ytilde = X_train @ beta\n",
"print(\"Training R2\")\n",
"print(R2(y_train,ytilde))\n",
"print(\"Training MSE\")\n",
"print(MSE(y_train,ytilde))\n",
"ypredict = X_test @ beta\n",
"print(\"Test R2\")\n",
"print(R2(y_test,ypredict))\n",
"print(\"Test MSE\")\n",
"print(MSE(y_test,ypredict))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"!esol\n",
"\n",
"\n",
"\n",
@@ -3545,9 +3485,7 @@
"\n",
"!bsubex\n",
"Add now a model which allows you to make polynomials up to degree $15$. Perform a standard OLS fitting of the training data and compute the MSE and $R2$ for the training and test data and plot both test and training data MSE and $R2$ as functions of the polynomial degree. Compare what you see with Figure 2.11 of Hastie et al. Comment your results. For which polynomial degree do you find an optimal MSE (smallest value)?\n",
"!bsol\n",
"Here you simply need to change the degree of the polynomial in the above code to $n=15$.\n",
"!esol\n",
"\n",
"!esubex"
]
}
+458 -73
View File
@@ -781,13 +781,13 @@ example of the functionality of <strong>Scikit-Learn</strong>.</p>
</div>
<div class="cell_output docutils container">
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>The intercept alpha:
[2.18954529]
[1.98452685]
Coefficient beta :
[[4.57051369]]
Mean squared error: 0.23
Variance score: 0.88
[[5.00109273]]
Mean squared error: 0.19
Variance score: 0.92
Mean squared log error: 0.01
Mean absolute error: 0.39
Mean absolute error: 0.35
</pre></div>
</div>
<img alt="_images/chapter1_13_1.png" src="_images/chapter1_13_1.png" />
@@ -888,7 +888,7 @@ a linear <span class="math notranslate nohighlight">\(x\)</span>-dependence we s
</div>
<div class="cell_output docutils container">
<img alt="_images/chapter1_27_0.png" src="_images/chapter1_27_0.png" />
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>0.004999999999999997
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>0.004999999999999996
</pre></div>
</div>
</div>
@@ -997,18 +997,6 @@ After having downloaded this file to our own computer, we are now ready to read
</pre></div>
</div>
</div>
<div class="cell_output docutils container">
<div class="output traceback highlight-ipythontb notranslate"><div class="highlight"><pre><span></span><span class="gt">---------------------------------------------------------------------------</span>
<span class="ne">FileNotFoundError</span><span class="g g-Whitespace"> </span>Traceback (most recent call last)
<span class="o">&lt;</span><span class="n">ipython</span><span class="o">-</span><span class="nb">input</span><span class="o">-</span><span class="mi">5</span><span class="o">-</span><span class="mi">3</span><span class="n">cd19a0768e1</span><span class="o">&gt;</span> <span class="ow">in</span> <span class="o">&lt;</span><span class="n">module</span><span class="o">&gt;</span>
<span class="g g-Whitespace"> </span><span class="mi">31</span> <span class="n">plt</span><span class="o">.</span><span class="n">savefig</span><span class="p">(</span><span class="n">image_path</span><span class="p">(</span><span class="n">fig_id</span><span class="p">)</span> <span class="o">+</span> <span class="s2">&quot;.png&quot;</span><span class="p">,</span> <span class="nb">format</span><span class="o">=</span><span class="s1">&#39;png&#39;</span><span class="p">)</span>
<span class="g g-Whitespace"> </span><span class="mi">32</span>
<span class="ne">---&gt; </span><span class="mi">33</span> <span class="n">infile</span> <span class="o">=</span> <span class="nb">open</span><span class="p">(</span><span class="n">data_path</span><span class="p">(</span><span class="s2">&quot;MassEval2016.dat&quot;</span><span class="p">),</span><span class="s1">&#39;r&#39;</span><span class="p">)</span>
<span class="ne">FileNotFoundError</span>: [Errno 2] No such file or directory: &#39;DataFiles/MassEval2016.dat&#39;
</pre></div>
</div>
</div>
</div>
<p>Before we proceed, we define also a function for making our plots. You can obviously avoid this and simply set up various <strong>matplotlib</strong> commands every time you need them. You may however find it convenient to collect all such commands in one function and simply call this function.</p>
<div class="cell docutils container">
@@ -1049,6 +1037,11 @@ data) to actually open the file and simply take a look at it!</p>
</pre></div>
</div>
</div>
<div class="cell_output docutils container">
<div class="output text_plain highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>&#39; \nThis is taken from the data file of the mass 2016 evaluation. \nAll files are 3436 lines long with 124 character per line. \n Headers are 39 lines long. \n col 1 : Fortran character control: 1 = page feed 0 = line feed \n format : a1,i3,i5,i5,i5,1x,a3,a4,1x,f13.5,f11.5,f11.3,f9.3,1x,a2,f11.3,f9.3,1x,i3,1x,f12.5,f11.5 \n These formats are reflected in the pandas widths variable below, see the statement \n widths=(1,3,5,5,5,1,3,4,1,13,11,11,9,1,2,11,9,1,3,1,12,11,1), \n Pandas has also a variable header, with length 39 in this case. \n&#39;
</pre></div>
</div>
</div>
</div>
<p>The data we are interested in are in columns 2, 3, 4 and 11, giving us
the number of neutrons, protons, mass numbers and binding energies,
@@ -1098,6 +1091,25 @@ the number of nucleons <span class="math notranslate nohighlight">\(A\)</span>,
</pre></div>
</div>
</div>
<div class="cell_output docutils container">
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span> N Z A Element Ebinding
A
1 0 0 1 1 H 0.000000
2 1 1 1 2 H 1.112283
3 2 2 1 3 H 2.827265
4 6 2 2 4 He 7.073915
5 9 3 2 5 He 5.512132
... ... ... ... ... ...
264 3304 156 108 264 Hs 7.298375
265 3310 157 108 265 Hs 7.296247
266 3317 158 108 266 Hs 7.298273
269 3338 159 110 269 Ds 7.250154
270 3344 160 110 270 Ds 7.253775
[267 rows x 5 columns]
</pre></div>
</div>
</div>
</div>
<p>The next step, and we will define this mathematically later, is to set up the so-called <strong>design matrix</strong>. We will throughout call this matrix <span class="math notranslate nohighlight">\(\boldsymbol{X}\)</span>.
It has dimensionality <span class="math notranslate nohighlight">\(p\times n\)</span>, where <span class="math notranslate nohighlight">\(n\)</span> is the number of data points and <span class="math notranslate nohighlight">\(p\)</span> are the so-called predictors. In our case here they are given by the number of polynomials in <span class="math notranslate nohighlight">\(A\)</span> we wish to include in the fit.</p>
@@ -1150,6 +1162,16 @@ Now we can print measures of how our fit is doing, the coefficients from the fit
</pre></div>
</div>
</div>
<div class="cell_output docutils container">
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>Mean squared error: 0.04
Variance score: 0.95
Mean absolute error: 0.05
[ 0.00000000e+00 7.06492086e-03 -1.73091052e-01 -1.66020213e+01
1.17385778e+00] 15.212327334149492
</pre></div>
</div>
<img alt="_images/chapter1_57_1.png" src="_images/chapter1_57_1.png" />
</div>
</div>
<p>As a teaser, let us now see how we can do this with decision trees using <strong>scikit-learn</strong>. Later we will switch to so-called <strong>random forests</strong>!</p>
<div class="cell docutils container">
@@ -1186,6 +1208,27 @@ Now we can print measures of how our fit is doing, the coefficients from the fit
</pre></div>
</div>
</div>
<div class="cell_output docutils container">
<img alt="_images/chapter1_59_0.png" src="_images/chapter1_59_0.png" />
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span> N Z A Element Ebinding Eapprox
A
1 0 0 1 1 H 0.000000 0.000000
2 1 1 1 2 H 1.112283 1.112283
3 2 2 1 3 H 2.827265 2.827265
4 6 2 2 4 He 7.073915 7.073915
5 9 3 2 5 He 5.512132 5.512132
... ... ... ... ... ... ...
264 3304 156 108 264 Hs 7.298375 7.298375
265 3310 157 108 265 Hs 7.296247 7.297260
266 3317 158 108 266 Hs 7.298273 7.297260
269 3338 159 110 269 Ds 7.250154 7.250154
270 3344 160 110 270 Ds 7.253775 7.253775
[267 rows x 6 columns]
0.009883615646716184
</pre></div>
</div>
</div>
</div>
<p>The <strong>seaborn</strong> package allows us to visualize data in an efficient way. Note that we use <strong>scikit-learn</strong>s multi-layer perceptron (or feed forward neural network)
functionality.</p>
@@ -1223,6 +1266,91 @@ functionality.</p>
</pre></div>
</div>
</div>
<div class="cell_output docutils container">
<div class="output stderr highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/sklearn/neural_network/_multilayer_perceptron.py:582: ConvergenceWarning: Stochastic Optimizer: Maximum iterations (100) reached and the optimization hasn&#39;t converged yet.
warnings.warn(
/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/sklearn/neural_network/_multilayer_perceptron.py:582: ConvergenceWarning: Stochastic Optimizer: Maximum iterations (100) reached and the optimization hasn&#39;t converged yet.
warnings.warn(
/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/sklearn/neural_network/_multilayer_perceptron.py:582: ConvergenceWarning: Stochastic Optimizer: Maximum iterations (100) reached and the optimization hasn&#39;t converged yet.
warnings.warn(
/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/sklearn/neural_network/_multilayer_perceptron.py:582: ConvergenceWarning: Stochastic Optimizer: Maximum iterations (100) reached and the optimization hasn&#39;t converged yet.
warnings.warn(
</pre></div>
</div>
<div class="output stderr highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/sklearn/neural_network/_multilayer_perceptron.py:582: ConvergenceWarning: Stochastic Optimizer: Maximum iterations (100) reached and the optimization hasn&#39;t converged yet.
warnings.warn(
/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/sklearn/neural_network/_multilayer_perceptron.py:582: ConvergenceWarning: Stochastic Optimizer: Maximum iterations (100) reached and the optimization hasn&#39;t converged yet.
warnings.warn(
/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/sklearn/neural_network/_multilayer_perceptron.py:582: ConvergenceWarning: Stochastic Optimizer: Maximum iterations (100) reached and the optimization hasn&#39;t converged yet.
warnings.warn(
/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/sklearn/neural_network/_multilayer_perceptron.py:582: ConvergenceWarning: Stochastic Optimizer: Maximum iterations (100) reached and the optimization hasn&#39;t converged yet.
warnings.warn(
</pre></div>
</div>
<div class="output stderr highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/sklearn/neural_network/_multilayer_perceptron.py:582: ConvergenceWarning: Stochastic Optimizer: Maximum iterations (100) reached and the optimization hasn&#39;t converged yet.
warnings.warn(
/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/sklearn/neural_network/_multilayer_perceptron.py:582: ConvergenceWarning: Stochastic Optimizer: Maximum iterations (100) reached and the optimization hasn&#39;t converged yet.
warnings.warn(
/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/sklearn/neural_network/_multilayer_perceptron.py:582: ConvergenceWarning: Stochastic Optimizer: Maximum iterations (100) reached and the optimization hasn&#39;t converged yet.
warnings.warn(
/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/sklearn/neural_network/_multilayer_perceptron.py:582: ConvergenceWarning: Stochastic Optimizer: Maximum iterations (100) reached and the optimization hasn&#39;t converged yet.
warnings.warn(
</pre></div>
</div>
<div class="output stderr highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/sklearn/neural_network/_multilayer_perceptron.py:582: ConvergenceWarning: Stochastic Optimizer: Maximum iterations (100) reached and the optimization hasn&#39;t converged yet.
warnings.warn(
/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/sklearn/neural_network/_multilayer_perceptron.py:582: ConvergenceWarning: Stochastic Optimizer: Maximum iterations (100) reached and the optimization hasn&#39;t converged yet.
warnings.warn(
/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/sklearn/neural_network/_multilayer_perceptron.py:582: ConvergenceWarning: Stochastic Optimizer: Maximum iterations (100) reached and the optimization hasn&#39;t converged yet.
warnings.warn(
/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/sklearn/neural_network/_multilayer_perceptron.py:582: ConvergenceWarning: Stochastic Optimizer: Maximum iterations (100) reached and the optimization hasn&#39;t converged yet.
warnings.warn(
</pre></div>
</div>
<div class="output stderr highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/sklearn/neural_network/_multilayer_perceptron.py:582: ConvergenceWarning: Stochastic Optimizer: Maximum iterations (100) reached and the optimization hasn&#39;t converged yet.
warnings.warn(
/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/sklearn/neural_network/_multilayer_perceptron.py:582: ConvergenceWarning: Stochastic Optimizer: Maximum iterations (100) reached and the optimization hasn&#39;t converged yet.
warnings.warn(
/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/sklearn/neural_network/_multilayer_perceptron.py:582: ConvergenceWarning: Stochastic Optimizer: Maximum iterations (100) reached and the optimization hasn&#39;t converged yet.
warnings.warn(
/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/sklearn/neural_network/_multilayer_perceptron.py:582: ConvergenceWarning: Stochastic Optimizer: Maximum iterations (100) reached and the optimization hasn&#39;t converged yet.
warnings.warn(
</pre></div>
</div>
<div class="output stderr highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/sklearn/neural_network/_multilayer_perceptron.py:582: ConvergenceWarning: Stochastic Optimizer: Maximum iterations (100) reached and the optimization hasn&#39;t converged yet.
warnings.warn(
/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/sklearn/neural_network/_multilayer_perceptron.py:582: ConvergenceWarning: Stochastic Optimizer: Maximum iterations (100) reached and the optimization hasn&#39;t converged yet.
warnings.warn(
/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/sklearn/neural_network/_multilayer_perceptron.py:582: ConvergenceWarning: Stochastic Optimizer: Maximum iterations (100) reached and the optimization hasn&#39;t converged yet.
warnings.warn(
</pre></div>
</div>
<div class="output stderr highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/sklearn/neural_network/_multilayer_perceptron.py:582: ConvergenceWarning: Stochastic Optimizer: Maximum iterations (100) reached and the optimization hasn&#39;t converged yet.
warnings.warn(
/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/sklearn/neural_network/_multilayer_perceptron.py:582: ConvergenceWarning: Stochastic Optimizer: Maximum iterations (100) reached and the optimization hasn&#39;t converged yet.
warnings.warn(
/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/sklearn/neural_network/_multilayer_perceptron.py:582: ConvergenceWarning: Stochastic Optimizer: Maximum iterations (100) reached and the optimization hasn&#39;t converged yet.
warnings.warn(
/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/sklearn/neural_network/_multilayer_perceptron.py:582: ConvergenceWarning: Stochastic Optimizer: Maximum iterations (100) reached and the optimization hasn&#39;t converged yet.
warnings.warn(
</pre></div>
</div>
<div class="output stderr highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/sklearn/neural_network/_multilayer_perceptron.py:582: ConvergenceWarning: Stochastic Optimizer: Maximum iterations (100) reached and the optimization hasn&#39;t converged yet.
warnings.warn(
/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/sklearn/neural_network/_multilayer_perceptron.py:582: ConvergenceWarning: Stochastic Optimizer: Maximum iterations (100) reached and the optimization hasn&#39;t converged yet.
warnings.warn(
</pre></div>
</div>
<div class="output stderr highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/sklearn/neural_network/_multilayer_perceptron.py:582: ConvergenceWarning: Stochastic Optimizer: Maximum iterations (100) reached and the optimization hasn&#39;t converged yet.
warnings.warn(
/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/sklearn/neural_network/_multilayer_perceptron.py:582: ConvergenceWarning: Stochastic Optimizer: Maximum iterations (100) reached and the optimization hasn&#39;t converged yet.
warnings.warn(
/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/sklearn/neural_network/_multilayer_perceptron.py:582: ConvergenceWarning: Stochastic Optimizer: Maximum iterations (100) reached and the optimization hasn&#39;t converged yet.
warnings.warn(
</pre></div>
</div>
<img alt="_images/chapter1_61_9.png" src="_images/chapter1_61_9.png" />
</div>
</div>
</div>
</div>
@@ -1459,6 +1587,133 @@ our matrix as <span class="math notranslate nohighlight">\(\boldsymbol{X}\in {\m
</pre></div>
</div>
</div>
<div class="cell_output docutils container">
<div class="output text_html"><div>
<style scoped>
.dataframe tbody tr th:only-of-type {
vertical-align: middle;
}
.dataframe tbody tr th {
vertical-align: top;
}
.dataframe thead th {
text-align: right;
}
</style>
<table border="1" class="dataframe">
<thead>
<tr style="text-align: right;">
<th></th>
<th>1</th>
<th>A</th>
<th>A^(2/3)</th>
<th>A^(-1/3)</th>
<th>1/A</th>
</tr>
<tr>
<th>A</th>
<th></th>
<th></th>
<th></th>
<th></th>
<th></th>
</tr>
</thead>
<tbody>
<tr>
<th>1</th>
<td>1.0</td>
<td>1.0</td>
<td>1.000000</td>
<td>1.000000</td>
<td>1.000000</td>
</tr>
<tr>
<th>2</th>
<td>1.0</td>
<td>2.0</td>
<td>1.587401</td>
<td>0.793701</td>
<td>0.500000</td>
</tr>
<tr>
<th>3</th>
<td>1.0</td>
<td>3.0</td>
<td>2.080084</td>
<td>0.693361</td>
<td>0.333333</td>
</tr>
<tr>
<th>4</th>
<td>1.0</td>
<td>4.0</td>
<td>2.519842</td>
<td>0.629961</td>
<td>0.250000</td>
</tr>
<tr>
<th>5</th>
<td>1.0</td>
<td>5.0</td>
<td>2.924018</td>
<td>0.584804</td>
<td>0.200000</td>
</tr>
<tr>
<th>...</th>
<td>...</td>
<td>...</td>
<td>...</td>
<td>...</td>
<td>...</td>
</tr>
<tr>
<th>264</th>
<td>1.0</td>
<td>264.0</td>
<td>41.153106</td>
<td>0.155883</td>
<td>0.003788</td>
</tr>
<tr>
<th>265</th>
<td>1.0</td>
<td>265.0</td>
<td>41.256962</td>
<td>0.155687</td>
<td>0.003774</td>
</tr>
<tr>
<th>266</th>
<td>1.0</td>
<td>266.0</td>
<td>41.360688</td>
<td>0.155491</td>
<td>0.003759</td>
</tr>
<tr>
<th>269</th>
<td>1.0</td>
<td>269.0</td>
<td>41.671089</td>
<td>0.154911</td>
<td>0.003717</td>
</tr>
<tr>
<th>270</th>
<td>1.0</td>
<td>270.0</td>
<td>41.774300</td>
<td>0.154720</td>
<td>0.003704</td>
</tr>
</tbody>
</table>
<p>267 rows × 5 columns</p>
</div></div></div>
</div>
<p>With <span class="math notranslate nohighlight">\(\boldsymbol{\beta}\in {\mathbb{R}}^{p\times 1}\)</span>, it means that we will hereafter write our equations for the approximation as</p>
<div class="math notranslate nohighlight">
@@ -1610,7 +1865,7 @@ C
K</p>
<div class="math notranslate nohighlight">
\[
\frac{\partial \log{\vert\boldsymbol{A}\vert}}{\partial \boldsymbol{A}} = (\boldsymbol{A}^{-1})^T.
\frac{\partial\log{\vert\boldsymbol{A}\vert}}{\partial \boldsymbol{A}}=(\boldsymbol{A}^{-1})^T.
\]</div>
<p>The residuals <span class="math notranslate nohighlight">\(\boldsymbol{\epsilon}\)</span> are in turn given by</p>
<div class="math notranslate nohighlight">
@@ -1668,6 +1923,9 @@ write</p>
</pre></div>
</div>
</div>
<div class="cell_output docutils container">
<img alt="_images/chapter1_132_0.png" src="_images/chapter1_132_0.png" />
</div>
</div>
<p>We can easily test our fit by computing the <span class="math notranslate nohighlight">\(R2\)</span> score that we discussed in connection with the functionality of <strong>Scikit-Learn</strong> in the introductory slides.
Since we are not using <strong>Scikit-Learn</strong> here we can define our own <span class="math notranslate nohighlight">\(R2\)</span> function as</p>
@@ -1686,6 +1944,11 @@ Since we are not using <strong>Scikit-Learn</strong> here we can define our own
</pre></div>
</div>
</div>
<div class="cell_output docutils container">
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>0.9547578478889096
</pre></div>
</div>
</div>
</div>
<p>We can easily add our <strong>MSE</strong> score as</p>
<div class="cell docutils container">
@@ -1698,6 +1961,11 @@ Since we are not using <strong>Scikit-Learn</strong> here we can define our own
</pre></div>
</div>
</div>
<div class="cell_output docutils container">
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>0.03787596148305236
</pre></div>
</div>
</div>
</div>
<p>and finally the relative error as</p>
<div class="cell docutils container">
@@ -1708,6 +1976,23 @@ Since we are not using <strong>Scikit-Learn</strong> here we can define our own
</pre></div>
</div>
</div>
<div class="cell_output docutils container">
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>A
1 0 inf
2 1 1.123190
3 2 0.327631
4 6 0.344172
5 9 0.044402
...
264 3304 0.009911
265 3310 0.009154
266 3317 0.007824
269 3338 0.011347
270 3344 0.009790
Name: Ebinding, Length: 267, dtype: float64
</pre></div>
</div>
</div>
</div>
<div class="section" id="the-chi-2-function">
<h3><span class="section-number">3.4.1. </span>The <span class="math notranslate nohighlight">\(\chi^2\)</span> function<a class="headerlink" href="#the-chi-2-function" title="Permalink to this headline"></a></h3>
@@ -1935,6 +2220,19 @@ hyperparameter <span class="math notranslate nohighlight">\(\lambda\)</span>, al
</pre></div>
</div>
</div>
<div class="cell_output docutils container">
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>Mean squared error: 12.36
Variance score: 1.00
Mean absolute error: 2.83
[ 0. 618.32047562 -861.13519106 1404.91549644] -11.057088709963637
Mean squared error: 197.93
Variance score: 1.00
Mean absolute error: 11.69
[ 0. 28.18220995 282.79902342 842.30879705] 12.946893955211749
</pre></div>
</div>
<img alt="_images/chapter1_179_1.png" src="_images/chapter1_179_1.png" />
</div>
</div>
<p>The above simple polynomial in density <span class="math notranslate nohighlight">\(\rho\)</span> gives an excellent fit
to the data.</p>
@@ -2024,6 +2322,18 @@ but now splitting the data into a training set and a test set.</p>
</pre></div>
</div>
</div>
<div class="cell_output docutils container">
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>Training R2
0.9999864543345858
Training MSE
6.180092462880674
Test R2
0.9999822527140678
Test MSE
7.205466494327873
</pre></div>
</div>
</div>
</div>
</div>
<div class="section" id="the-boston-housing-data-example">
@@ -2078,6 +2388,11 @@ the house using the features (predictors) listed here.</p>
</pre></div>
</div>
</div>
<div class="cell_output docutils container">
<div class="output text_plain highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>dict_keys([&#39;data&#39;, &#39;target&#39;, &#39;feature_names&#39;, &#39;DESCR&#39;, &#39;filename&#39;])
</pre></div>
</div>
</div>
</div>
<p>Then we invoke Pandas</p>
<div class="cell docutils container">
@@ -2097,6 +2412,25 @@ the house using the features (predictors) listed here.</p>
</pre></div>
</div>
</div>
<div class="cell_output docutils container">
<div class="output text_plain highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>CRIM 0
ZN 0
INDUS 0
CHAS 0
NOX 0
RM 0
AGE 0
DIS 0
RAD 0
TAX 0
PTRATIO 0
B 0
LSTAT 0
MEDV 0
dtype: int64
</pre></div>
</div>
</div>
</div>
<p>We can then visualize the data</p>
<div class="cell docutils container">
@@ -2110,6 +2444,13 @@ the house using the features (predictors) listed here.</p>
</pre></div>
</div>
</div>
<div class="cell_output docutils container">
<div class="output stderr highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/seaborn/distributions.py:2551: FutureWarning: `distplot` is a deprecated function and will be removed in a future version. Please adapt your code to use either `displot` (a figure-level function with similar flexibility) or `histplot` (an axes-level function for histograms).
warnings.warn(msg, FutureWarning)
</pre></div>
</div>
<img alt="_images/chapter1_191_1.png" src="_images/chapter1_191_1.png" />
</div>
</div>
<p>It is now useful to look at the correlation matrix</p>
<div class="cell docutils container">
@@ -2122,6 +2463,12 @@ the house using the features (predictors) listed here.</p>
</pre></div>
</div>
</div>
<div class="cell_output docutils container">
<div class="output text_plain highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>&lt;AxesSubplot:&gt;
</pre></div>
</div>
<img alt="_images/chapter1_193_1.png" src="_images/chapter1_193_1.png" />
</div>
</div>
<p>From the above coorelation plot we can see that <strong>MEDV</strong> is strongly correlated to <strong>LSTAT</strong> and <strong>RM</strong>. We see also that <strong>RAD</strong> and <strong>TAX</strong> are stronly correlated, but we dont include this in our features together to avoid multi-colinearity</p>
<div class="cell docutils container">
@@ -2142,6 +2489,9 @@ the house using the features (predictors) listed here.</p>
</pre></div>
</div>
</div>
<div class="cell_output docutils container">
<img alt="_images/chapter1_195_0.png" src="_images/chapter1_195_0.png" />
</div>
</div>
<p>Now we start training our model</p>
<div class="cell docutils container">
@@ -2167,6 +2517,14 @@ the house using the features (predictors) listed here.</p>
</pre></div>
</div>
</div>
<div class="cell_output docutils container">
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>(404, 2)
(102, 2)
(404,)
(102,)
</pre></div>
</div>
</div>
</div>
<p>Then we use the linear regression functionality from <strong>Scikit-Learn</strong></p>
<div class="cell docutils container">
@@ -2205,6 +2563,20 @@ the house using the features (predictors) listed here.</p>
</pre></div>
</div>
</div>
<div class="cell_output docutils container">
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>The model performance for training set
--------------------------------------
RMSE is 5.637129335071195
R2 score is 0.6300745149331701
The model performance for testing set
--------------------------------------
RMSE is 5.137400784702912
R2 score is 0.6628996975186952
</pre></div>
</div>
</div>
</div>
<div class="cell docutils container">
<div class="cell_input docutils container">
@@ -2215,6 +2587,9 @@ the house using the features (predictors) listed here.</p>
</pre></div>
</div>
</div>
<div class="cell_output docutils container">
<img alt="_images/chapter1_202_0.png" src="_images/chapter1_202_0.png" />
</div>
</div>
</div>
<div class="section" id="reducing-the-number-of-degrees-of-freedom-overarching-view">
@@ -2364,6 +2739,36 @@ techniques.</p>
</pre></div>
</div>
</div>
<div class="cell_output docutils container">
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>MSE before scaling: 0.00
R2 score before scaling 0.99
Feature min values before scaling:
[1.00000000e+00 1.10094646e-03 9.51523276e-04 1.21208310e-06
1.04757618e-06 9.05396545e-07 1.33443859e-09 1.15332528e-09
9.96793117e-10 8.61505887e-10 1.46914544e-12 1.26974938e-12
1.09741585e-12 9.48471852e-13 8.19742904e-13 1.61745046e-15
1.39792608e-15 1.20819609e-15 1.04421672e-15 9.02493044e-16
7.80004454e-16]
Feature max values before scaling:
[1. 0.99825997 0.99883879 0.99652296 0.99710078 0.99767893
0.99478898 0.9953658 0.99594294 0.99652042 0.99305802 0.99363383
0.99420997 0.99478645 0.99536326 0.99133007 0.99190487 0.99248001
0.99305549 0.99363129 0.99420743]
Feature min values after scaling:
[ 0. -1.61869821 -1.66880047 -1.06209126 -1.08170444 -1.10131725
-0.84087101 -0.85396354 -0.86692943 -0.87972591 -0.71351486 -0.7240496
-0.73453972 -0.74495014 -0.75524378 -0.62783293 -0.63680118 -0.64580686
-0.65482578 -0.66383151 -0.67279536]
Feature max values after scaling:
[0. 1.78944806 1.68342382 2.34172919 2.25457052 2.16553696
2.7899453 2.71281409 2.63374631 2.55280484 3.17117385 3.10307631
3.03303359 2.96104648 2.88712946 3.50394742 3.44395541 3.38216436
3.31853484 3.25303483 3.1856411 ]
MSE after scaling: 0.00
R2 score for scaled data: 0.99
</pre></div>
</div>
</div>
</div>
</div>
</div>
@@ -2460,55 +2865,6 @@ R^2(\hat{y}, \tilde{\hat{y}}) = 1 - \frac{\sum_{i=0}^{n - 1} (y_i - \tilde{y}_i)
\]</div>
<p>You can use the functionality included in scikit-learn. If you feel for it, you can use your own program and define functions which compute the above two functions.
Discuss the meaning of these results. Try also to vary the coefficient in front of the added stochastic noise term and discuss the quality of the fits.</p>
<p>!bsol
The code here is an example of where we define our own design matrix and fit parameters <span class="math notranslate nohighlight">\(\beta\)</span>.</p>
<div class="cell docutils container">
<div class="cell_input docutils container">
<div class="highlight-ipython3 notranslate"><div class="highlight"><pre><span></span><span class="kn">import</span> <span class="nn">os</span>
<span class="kn">import</span> <span class="nn">numpy</span> <span class="k">as</span> <span class="nn">np</span>
<span class="kn">import</span> <span class="nn">pandas</span> <span class="k">as</span> <span class="nn">pd</span>
<span class="kn">import</span> <span class="nn">matplotlib.pyplot</span> <span class="k">as</span> <span class="nn">plt</span>
<span class="kn">from</span> <span class="nn">sklearn.model_selection</span> <span class="kn">import</span> <span class="n">train_test_split</span>
<span class="k">def</span> <span class="nf">save_fig</span><span class="p">(</span><span class="n">fig_id</span><span class="p">):</span>
<span class="n">plt</span><span class="o">.</span><span class="n">savefig</span><span class="p">(</span><span class="n">image_path</span><span class="p">(</span><span class="n">fig_id</span><span class="p">)</span> <span class="o">+</span> <span class="s2">&quot;.png&quot;</span><span class="p">,</span> <span class="nb">format</span><span class="o">=</span><span class="s1">&#39;png&#39;</span><span class="p">)</span>
<span class="k">def</span> <span class="nf">R2</span><span class="p">(</span><span class="n">y_data</span><span class="p">,</span> <span class="n">y_model</span><span class="p">):</span>
<span class="k">return</span> <span class="mi">1</span> <span class="o">-</span> <span class="n">np</span><span class="o">.</span><span class="n">sum</span><span class="p">((</span><span class="n">y_data</span> <span class="o">-</span> <span class="n">y_model</span><span class="p">)</span> <span class="o">**</span> <span class="mi">2</span><span class="p">)</span> <span class="o">/</span> <span class="n">np</span><span class="o">.</span><span class="n">sum</span><span class="p">((</span><span class="n">y_data</span> <span class="o">-</span> <span class="n">np</span><span class="o">.</span><span class="n">mean</span><span class="p">(</span><span class="n">y_data</span><span class="p">))</span> <span class="o">**</span> <span class="mi">2</span><span class="p">)</span>
<span class="k">def</span> <span class="nf">MSE</span><span class="p">(</span><span class="n">y_data</span><span class="p">,</span><span class="n">y_model</span><span class="p">):</span>
<span class="n">n</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">size</span><span class="p">(</span><span class="n">y_model</span><span class="p">)</span>
<span class="k">return</span> <span class="n">np</span><span class="o">.</span><span class="n">sum</span><span class="p">((</span><span class="n">y_data</span><span class="o">-</span><span class="n">y_model</span><span class="p">)</span><span class="o">**</span><span class="mi">2</span><span class="p">)</span><span class="o">/</span><span class="n">n</span>
<span class="n">x</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">random</span><span class="o">.</span><span class="n">rand</span><span class="p">(</span><span class="mi">100</span><span class="p">)</span>
<span class="n">y</span> <span class="o">=</span> <span class="mf">2.0</span><span class="o">+</span><span class="mi">5</span><span class="o">*</span><span class="n">x</span><span class="o">*</span><span class="n">x</span><span class="o">+</span><span class="mf">0.1</span><span class="o">*</span><span class="n">np</span><span class="o">.</span><span class="n">random</span><span class="o">.</span><span class="n">randn</span><span class="p">(</span><span class="mi">100</span><span class="p">)</span>
<span class="c1"># The design matrix now as function of a given polynomial</span>
<span class="n">X</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">((</span><span class="nb">len</span><span class="p">(</span><span class="n">x</span><span class="p">),</span><span class="mi">3</span><span class="p">))</span>
<span class="n">X</span><span class="p">[:,</span><span class="mi">0</span><span class="p">]</span> <span class="o">=</span> <span class="mf">1.0</span>
<span class="n">X</span><span class="p">[:,</span><span class="mi">1</span><span class="p">]</span> <span class="o">=</span> <span class="n">x</span>
<span class="n">X</span><span class="p">[:,</span><span class="mi">2</span><span class="p">]</span> <span class="o">=</span> <span class="n">x</span><span class="o">**</span><span class="mi">2</span>
<span class="c1"># We split the data in test and training data</span>
<span class="n">X_train</span><span class="p">,</span> <span class="n">X_test</span><span class="p">,</span> <span class="n">y_train</span><span class="p">,</span> <span class="n">y_test</span> <span class="o">=</span> <span class="n">train_test_split</span><span class="p">(</span><span class="n">X</span><span class="p">,</span> <span class="n">y</span><span class="p">,</span> <span class="n">test_size</span><span class="o">=</span><span class="mf">0.2</span><span class="p">)</span>
<span class="c1"># matrix inversion to find beta</span>
<span class="n">beta</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">linalg</span><span class="o">.</span><span class="n">inv</span><span class="p">(</span><span class="n">X_train</span><span class="o">.</span><span class="n">T</span> <span class="o">@</span> <span class="n">X_train</span><span class="p">)</span> <span class="o">@</span> <span class="n">X_train</span><span class="o">.</span><span class="n">T</span> <span class="o">@</span> <span class="n">y_train</span>
<span class="nb">print</span><span class="p">(</span><span class="n">beta</span><span class="p">)</span>
<span class="c1"># and then make the prediction</span>
<span class="n">ytilde</span> <span class="o">=</span> <span class="n">X_train</span> <span class="o">@</span> <span class="n">beta</span>
<span class="nb">print</span><span class="p">(</span><span class="s2">&quot;Training R2&quot;</span><span class="p">)</span>
<span class="nb">print</span><span class="p">(</span><span class="n">R2</span><span class="p">(</span><span class="n">y_train</span><span class="p">,</span><span class="n">ytilde</span><span class="p">))</span>
<span class="nb">print</span><span class="p">(</span><span class="s2">&quot;Training MSE&quot;</span><span class="p">)</span>
<span class="nb">print</span><span class="p">(</span><span class="n">MSE</span><span class="p">(</span><span class="n">y_train</span><span class="p">,</span><span class="n">ytilde</span><span class="p">))</span>
<span class="n">ypredict</span> <span class="o">=</span> <span class="n">X_test</span> <span class="o">@</span> <span class="n">beta</span>
<span class="nb">print</span><span class="p">(</span><span class="s2">&quot;Test R2&quot;</span><span class="p">)</span>
<span class="nb">print</span><span class="p">(</span><span class="n">R2</span><span class="p">(</span><span class="n">y_test</span><span class="p">,</span><span class="n">ypredict</span><span class="p">))</span>
<span class="nb">print</span><span class="p">(</span><span class="s2">&quot;Test MSE&quot;</span><span class="p">)</span>
<span class="nb">print</span><span class="p">(</span><span class="n">MSE</span><span class="p">(</span><span class="n">y_test</span><span class="p">,</span><span class="n">ypredict</span><span class="p">))</span>
</pre></div>
</div>
</div>
</div>
<p>!esol</p>
</div>
<div class="section" id="exercise-normalizing-our-data">
<h3><span class="section-number">3.9.3. </span>Exercise: Normalizing our data<a class="headerlink" href="#exercise-normalizing-our-data" title="Permalink to this headline"></a></h3>
@@ -2549,6 +2905,38 @@ for testing. This can be done as follows with our design matrix <span class="mat
</pre></div>
</div>
</div>
<div class="cell_output docutils container">
<div class="output traceback highlight-ipythontb notranslate"><div class="highlight"><pre><span></span><span class="gt">---------------------------------------------------------------------------</span>
<span class="ne">ValueError</span><span class="g g-Whitespace"> </span>Traceback (most recent call last)
<span class="o">&lt;</span><span class="n">ipython</span><span class="o">-</span><span class="nb">input</span><span class="o">-</span><span class="mi">38</span><span class="o">-</span><span class="mi">9</span><span class="n">b9cf4fa1a95</span><span class="o">&gt;</span> <span class="ow">in</span> <span class="o">&lt;</span><span class="n">module</span><span class="o">&gt;</span>
<span class="g g-Whitespace"> </span><span class="mi">1</span> <span class="c1"># split in training and test data</span>
<span class="ne">----&gt; </span><span class="mi">2</span> <span class="n">X_train</span><span class="p">,</span> <span class="n">X_test</span><span class="p">,</span> <span class="n">y_train</span><span class="p">,</span> <span class="n">y_test</span> <span class="o">=</span> <span class="n">train_test_split</span><span class="p">(</span><span class="n">X</span><span class="p">,</span><span class="n">y</span><span class="p">,</span><span class="n">test_size</span><span class="o">=</span><span class="mf">0.2</span><span class="p">)</span>
<span class="nn">~/opt/anaconda3/lib/python3.8/site-packages/sklearn/model_selection/_split.py</span> in <span class="ni">train_test_split</span><span class="nt">(*arrays, **options)</span>
<span class="g g-Whitespace"> </span><span class="mi">2125</span> <span class="k">raise</span> <span class="ne">TypeError</span><span class="p">(</span><span class="s2">&quot;Invalid parameters passed: </span><span class="si">%s</span><span class="s2">&quot;</span> <span class="o">%</span> <span class="nb">str</span><span class="p">(</span><span class="n">options</span><span class="p">))</span>
<span class="g g-Whitespace"> </span><span class="mi">2126</span>
<span class="ne">-&gt; </span><span class="mi">2127</span> <span class="n">arrays</span> <span class="o">=</span> <span class="n">indexable</span><span class="p">(</span><span class="o">*</span><span class="n">arrays</span><span class="p">)</span>
<span class="g g-Whitespace"> </span><span class="mi">2128</span>
<span class="g g-Whitespace"> </span><span class="mi">2129</span> <span class="n">n_samples</span> <span class="o">=</span> <span class="n">_num_samples</span><span class="p">(</span><span class="n">arrays</span><span class="p">[</span><span class="mi">0</span><span class="p">])</span>
<span class="nn">~/opt/anaconda3/lib/python3.8/site-packages/sklearn/utils/validation.py</span> in <span class="ni">indexable</span><span class="nt">(*iterables)</span>
<span class="g g-Whitespace"> </span><span class="mi">290</span> <span class="s2">&quot;&quot;&quot;</span>
<span class="g g-Whitespace"> </span><span class="mi">291</span><span class="s2"> result = [_make_indexable(X) for X in iterables]</span>
<span class="ne">--&gt; </span><span class="mi">292</span><span class="s2"> check_consistent_length(*result)</span>
<span class="g g-Whitespace"> </span><span class="mi">293</span><span class="s2"> return result</span>
<span class="g g-Whitespace"> </span><span class="mi">294</span><span class="s2"> </span>
<span class="nn">~/opt/anaconda3/lib/python3.8/site-packages/sklearn/utils/validation.py</span> in <span class="ni">check_consistent_length</span><span class="nt">(*arrays)</span>
<span class="g g-Whitespace"> </span><span class="mi">253</span><span class="s2"> uniques = np.unique(lengths)</span>
<span class="g g-Whitespace"> </span><span class="mi">254</span><span class="s2"> if len(uniques) &gt; 1:</span>
<span class="ne">--&gt; </span><span class="mi">255</span><span class="s2"> raise ValueError(&quot;Found input variables with inconsistent numbers of&quot;</span>
<span class="g g-Whitespace"> </span><span class="mi">256</span><span class="s2"> &quot; samples: </span><span class="si">%r</span><span class="s2">&quot; % [int(l) for l in lengths])</span>
<span class="g g-Whitespace"> </span><span class="mi">257</span><span class="s2"> </span>
<span class="ne">ValueError</span>: Found input variables with inconsistent numbers of samples: [1000, 100]
</pre></div>
</div>
</div>
</div>
<p>Then we can use the standard scaler to scale our data as</p>
<div class="cell docutils container">
@@ -2587,11 +2975,8 @@ Write a first code which sets up a design matrix <span class="math notranslate n
Perform an ordinary least squares and compute the means squared error and the <span class="math notranslate nohighlight">\(R2\)</span> factor for the training data and the test data, with and without scaling.
!esubex</p>
<p>!bsubex
Add now a model which allows you to make polynomials up to degree <span class="math notranslate nohighlight">\(15\)</span>. Perform a standard OLS fitting of the training data and compute the MSE and <span class="math notranslate nohighlight">\(R2\)</span> for the training and test data and plot both test and training data MSE and <span class="math notranslate nohighlight">\(R2\)</span> as functions of the polynomial degree. Compare what you see with Figure 2.11 of Hastie et al. Comment your results. For which polynomial degree do you find an optimal MSE (smallest value)?
!bsol
Here you simply need to change the degree of the polynomial in the above code to <span class="math notranslate nohighlight">\(n=15\)</span>.
!esol
!esubex</p>
Add now a model which allows you to make polynomials up to degree <span class="math notranslate nohighlight">\(15\)</span>. Perform a standard OLS fitting of the training data and compute the MSE and <span class="math notranslate nohighlight">\(R2\)</span> for the training and test data and plot both test and training data MSE and <span class="math notranslate nohighlight">\(R2\)</span> as functions of the polynomial degree. Compare what you see with Figure 2.11 of Hastie et al. Comment your results. For which polynomial degree do you find an optimal MSE (smallest value)?</p>
<p>!esubex</p>
</div>
</div>
</div>
@@ -0,0 +1,53 @@
Traceback (most recent call last):
File "/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/jupyter_cache/executors/utils.py", line 51, in single_nb_execution
executenb(
File "/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/nbclient/client.py", line 1087, in execute
return NotebookClient(nb=nb, resources=resources, km=km, **kwargs).execute()
File "/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/nbclient/util.py", line 74, in wrapped
return just_run(coro(*args, **kwargs))
File "/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/nbclient/util.py", line 53, in just_run
return loop.run_until_complete(coro)
File "/Users/mhjensen/opt/anaconda3/lib/python3.8/asyncio/base_events.py", line 616, in run_until_complete
return future.result()
File "/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/nbclient/client.py", line 540, in async_execute
await self.async_execute_cell(
File "/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/nbclient/client.py", line 832, in async_execute_cell
self._check_raise_for_error(cell, exec_reply)
File "/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/nbclient/client.py", line 740, in _check_raise_for_error
raise CellExecutionError.from_cell_and_msg(cell, exec_reply['content'])
nbclient.exceptions.CellExecutionError: An error occurred while executing the following cell:
------------------
# split in training and test data
X_train, X_test, y_train, y_test = train_test_split(X,y,test_size=0.2)
------------------
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
<ipython-input-38-9b9cf4fa1a95> in <module>
 1 # split in training and test data
----> 2 X_train, X_test, y_train, y_test = train_test_split(X,y,test_size=0.2)

~/opt/anaconda3/lib/python3.8/site-packages/sklearn/model_selection/_split.py in train_test_split(*arrays, **options)
 2125 raise TypeError("Invalid parameters passed: %s" % str(options))
 2126 
-> 2127 arrays = indexable(*arrays)
 2128 
 2129 n_samples = _num_samples(arrays[0])
~/opt/anaconda3/lib/python3.8/site-packages/sklearn/utils/validation.py in indexable(*iterables)
 290 """
 291 result = [_make_indexable(X) for X in iterables]
--> 292 check_consistent_length(*result)
 293 return result
 294 
~/opt/anaconda3/lib/python3.8/site-packages/sklearn/utils/validation.py in check_consistent_length(*arrays)
 253 uniques = np.unique(lengths)
 254 if len(uniques) > 1:
--> 255 raise ValueError("Found input variables with inconsistent numbers of"
 256 " samples: %r" % [int(l) for l in lengths])
 257 
ValueError: Found input variables with inconsistent numbers of samples: [1000, 100]
ValueError: Found input variables with inconsistent numbers of samples: [1000, 100]
@@ -0,0 +1,29 @@
Traceback (most recent call last):
File "/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/jupyter_cache/executors/utils.py", line 51, in single_nb_execution
executenb(
File "/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/nbclient/client.py", line 1087, in execute
return NotebookClient(nb=nb, resources=resources, km=km, **kwargs).execute()
File "/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/nbclient/util.py", line 74, in wrapped
return just_run(coro(*args, **kwargs))
File "/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/nbclient/util.py", line 53, in just_run
return loop.run_until_complete(coro)
File "/Users/mhjensen/opt/anaconda3/lib/python3.8/asyncio/base_events.py", line 616, in run_until_complete
return future.result()
File "/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/nbclient/client.py", line 540, in async_execute
await self.async_execute_cell(
File "/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/nbclient/client.py", line 832, in async_execute_cell
self._check_raise_for_error(cell, exec_reply)
File "/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/nbclient/client.py", line 740, in _check_raise_for_error
raise CellExecutionError.from_cell_and_msg(cell, exec_reply['content'])
nbclient.exceptions.CellExecutionError: An error occurred while executing the following cell:
------------------
pip3 install tensorflow
------------------
 File "<ipython-input-12-6ea927cc6e88>", line 1
 pip3 install tensorflow
 ^
SyntaxError: invalid syntax
SyntaxError: invalid syntax (<ipython-input-12-6ea927cc6e88>, line 1)
@@ -0,0 +1,104 @@
Traceback (most recent call last):
File "/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/jupyter_cache/executors/utils.py", line 51, in single_nb_execution
executenb(
File "/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/nbclient/client.py", line 1087, in execute
return NotebookClient(nb=nb, resources=resources, km=km, **kwargs).execute()
File "/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/nbclient/util.py", line 74, in wrapped
return just_run(coro(*args, **kwargs))
File "/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/nbclient/util.py", line 53, in just_run
return loop.run_until_complete(coro)
File "/Users/mhjensen/opt/anaconda3/lib/python3.8/asyncio/base_events.py", line 616, in run_until_complete
return future.result()
File "/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/nbclient/client.py", line 540, in async_execute
await self.async_execute_cell(
File "/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/nbclient/client.py", line 832, in async_execute_cell
self._check_raise_for_error(cell, exec_reply)
File "/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/nbclient/client.py", line 740, in _check_raise_for_error
raise CellExecutionError.from_cell_and_msg(cell, exec_reply['content'])
nbclient.exceptions.CellExecutionError: An error occurred while executing the following cell:
------------------
%matplotlib inline
from numpy import *
from numpy.random import randint, randn
from time import time
import matplotlib.mlab as mlab
import matplotlib.pyplot as plt
# Returns mean of bootstrap samples
def stat(data):
return mean(data)
# Bootstrap algorithm
def bootstrap(data, statistic, R):
t = zeros(R); n = len(data); inds = arange(n); t0 = time()
# non-parametric bootstrap
for i in range(R):
t[i] = statistic(data[randint(0,n,n)])
# analysis
print("Runtime: %g sec" % (time()-t0)); print("Bootstrap Statistics :")
print("original bias std. error")
print("%8g %8g %14g %15g" % (statistic(data), std(data),mean(t),std(t)))
return t
mu, sigma = 100, 15
datapoints = 10000
x = mu + sigma*random.randn(datapoints)
# bootstrap returns the data sample
t = bootstrap(x, stat, datapoints)
# the histogram of the bootstrapped data
n, binsboot, patches = plt.hist(t, 50, normed=1, facecolor='red', alpha=0.75)
# add a 'best fit' line
y = mlab.normpdf( binsboot, mean(t), std(t))
lt = plt.plot(binsboot, y, 'r--', linewidth=1)
plt.xlabel('Smarts')
plt.ylabel('Probability')
plt.axis([99.5, 100.6, 0, 3.0])
plt.grid(True)
plt.show()
------------------
---------------------------------------------------------------------------
AttributeError Traceback (most recent call last)
<ipython-input-2-772b904ae9cb> in <module>
 31 t = bootstrap(x, stat, datapoints)
 32 # the histogram of the bootstrapped data
---> 33 n, binsboot, patches = plt.hist(t, 50, normed=1, facecolor='red', alpha=0.75)
 34 
 35 # add a 'best fit' line
~/opt/anaconda3/lib/python3.8/site-packages/matplotlib/pyplot.py in hist(x, bins, range, density, weights, cumulative, bottom, histtype, align, orientation, rwidth, log, color, label, stacked, data, **kwargs)
 2683 orientation='vertical', rwidth=None, log=False, color=None,
 2684 label=None, stacked=False, *, data=None, **kwargs):
-> 2685 return gca().hist(
 2686 x, bins=bins, range=range, density=density, weights=weights,
 2687 cumulative=cumulative, bottom=bottom, histtype=histtype,
~/opt/anaconda3/lib/python3.8/site-packages/matplotlib/__init__.py in inner(ax, data, *args, **kwargs)
 1445 def inner(ax, *args, data=None, **kwargs):
 1446 if data is None:
-> 1447 return func(ax, *map(sanitize_sequence, args), **kwargs)
 1448 
 1449 bound = new_sig.bind(ax, *args, **kwargs)
~/opt/anaconda3/lib/python3.8/site-packages/matplotlib/axes/_axes.py in hist(self, x, bins, range, density, weights, cumulative, bottom, histtype, align, orientation, rwidth, log, color, label, stacked, **kwargs)
 6813 if patch:
 6814 p = patch[0]
-> 6815 p.update(kwargs)
 6816 if lbl is not None:
 6817 p.set_label(lbl)
~/opt/anaconda3/lib/python3.8/site-packages/matplotlib/artist.py in update(self, props)
 994 func = getattr(self, f"set_{k}", None)
 995 if not callable(func):
--> 996 raise AttributeError(f"{type(self).__name__!r} object "
 997 f"has no property {k!r}")
 998 ret.append(func(v))
AttributeError: 'Rectangle' object has no property 'normed'
AttributeError: 'Rectangle' object has no property 'normed'
@@ -0,0 +1,89 @@
Traceback (most recent call last):
File "/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/jupyter_cache/executors/utils.py", line 51, in single_nb_execution
executenb(
File "/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/nbclient/client.py", line 1087, in execute
return NotebookClient(nb=nb, resources=resources, km=km, **kwargs).execute()
File "/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/nbclient/util.py", line 74, in wrapped
return just_run(coro(*args, **kwargs))
File "/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/nbclient/util.py", line 53, in just_run
return loop.run_until_complete(coro)
File "/Users/mhjensen/opt/anaconda3/lib/python3.8/asyncio/base_events.py", line 616, in run_until_complete
return future.result()
File "/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/nbclient/client.py", line 540, in async_execute
await self.async_execute_cell(
File "/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/nbclient/client.py", line 832, in async_execute_cell
self._check_raise_for_error(cell, exec_reply)
File "/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/nbclient/client.py", line 740, in _check_raise_for_error
raise CellExecutionError.from_cell_and_msg(cell, exec_reply['content'])
nbclient.exceptions.CellExecutionError: An error occurred while executing the following cell:
------------------
%matplotlib inline
# Common imports
import os
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression, Ridge, Lasso
from sklearn.model_selection import train_test_split
from sklearn.utils import resample
from sklearn.metrics import mean_squared_error
from IPython.display import display
from pylab import plt, mpl
plt.style.use('seaborn')
mpl.rcParams['font.family'] = 'serif'
# Where to save the figures and data files
PROJECT_ROOT_DIR = "Results"
FIGURE_ID = "Results/FigureFiles"
DATA_ID = "DataFiles/"
if not os.path.exists(PROJECT_ROOT_DIR):
os.mkdir(PROJECT_ROOT_DIR)
if not os.path.exists(FIGURE_ID):
os.makedirs(FIGURE_ID)
if not os.path.exists(DATA_ID):
os.makedirs(DATA_ID)
def image_path(fig_id):
return os.path.join(FIGURE_ID, fig_id)
def data_path(dat_id):
return os.path.join(DATA_ID, dat_id)
def save_fig(fig_id):
plt.savefig(image_path(fig_id) + ".png", format='png')
infile = open(data_path("chddata.csv"),'r')
# Read the chd data as csv file and organize the data into arrays with age group, age, and chd
chd = pd.read_csv(infile, names=('ID', 'Age', 'Agegroup', 'CHD'))
chd.columns = ['ID', 'Age', 'Agegroup', 'CHD']
output = chd['CHD']
age = chd['Age']
agegroup = chd['Agegroup']
numberID = chd['ID']
display(chd)
plt.scatter(age, output, marker='o')
plt.axis([18,70.0,-0.1, 1.2])
plt.xlabel(r'Age')
plt.ylabel(r'CHD')
plt.title(r'Age distribution and Coronary heart disease')
plt.show()
------------------
---------------------------------------------------------------------------
FileNotFoundError Traceback (most recent call last)
<ipython-input-1-a77d5ac269b2> in <module>
 38 plt.savefig(image_path(fig_id) + ".png", format='png')
 39 
---> 40 infile = open(data_path("chddata.csv"),'r')
 41 
 42 # Read the chd data as csv file and organize the data into arrays with age group, age, and chd
FileNotFoundError: [Errno 2] No such file or directory: 'DataFiles/chddata.csv'
FileNotFoundError: [Errno 2] No such file or directory: 'DataFiles/chddata.csv'
@@ -0,0 +1,41 @@
Traceback (most recent call last):
File "/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/jupyter_cache/executors/utils.py", line 51, in single_nb_execution
executenb(
File "/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/nbclient/client.py", line 1087, in execute
return NotebookClient(nb=nb, resources=resources, km=km, **kwargs).execute()
File "/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/nbclient/util.py", line 74, in wrapped
return just_run(coro(*args, **kwargs))
File "/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/nbclient/util.py", line 53, in just_run
return loop.run_until_complete(coro)
File "/Users/mhjensen/opt/anaconda3/lib/python3.8/asyncio/base_events.py", line 616, in run_until_complete
return future.result()
File "/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/nbclient/client.py", line 540, in async_execute
await self.async_execute_cell(
File "/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/nbclient/client.py", line 832, in async_execute_cell
self._check_raise_for_error(cell, exec_reply)
File "/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/nbclient/client.py", line 740, in _check_raise_for_error
raise CellExecutionError.from_cell_and_msg(cell, exec_reply['content'])
nbclient.exceptions.CellExecutionError: An error occurred while executing the following cell:
------------------
# Import the necessary packages
import numpy
from cvxopt import matrix
from cvxopt import solvers
P = matrix(numpy.diag([1,0]), tc=d)
q = matrix(numpy.array([3,4]), tc=d)
G = matrix(numpy.array([[-1,0],[0,-1],[-1,-3],[2,5],[3,4]]), tc=d)
h = matrix(numpy.array([0,0,-15,100,80]), tc=d)
# Construct the QP, invoke solver
sol = solvers.qp(P,q,G,h)
# Extract optimal value and solution
sol[x]
sol[primal objective]
------------------
 File "<ipython-input-5-c46dd114b2af>", line 5
 P = matrix(numpy.diag([1,0]), tc=d)
 ^
SyntaxError: invalid character in identifier
SyntaxError: invalid character in identifier (<ipython-input-5-c46dd114b2af>, line 5)
@@ -0,0 +1,66 @@
Traceback (most recent call last):
File "/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/jupyter_cache/executors/utils.py", line 51, in single_nb_execution
executenb(
File "/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/nbclient/client.py", line 1087, in execute
return NotebookClient(nb=nb, resources=resources, km=km, **kwargs).execute()
File "/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/nbclient/util.py", line 74, in wrapped
return just_run(coro(*args, **kwargs))
File "/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/nbclient/util.py", line 53, in just_run
return loop.run_until_complete(coro)
File "/Users/mhjensen/opt/anaconda3/lib/python3.8/asyncio/base_events.py", line 616, in run_until_complete
return future.result()
File "/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/nbclient/client.py", line 540, in async_execute
await self.async_execute_cell(
File "/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/nbclient/client.py", line 832, in async_execute_cell
self._check_raise_for_error(cell, exec_reply)
File "/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/nbclient/client.py", line 740, in _check_raise_for_error
raise CellExecutionError.from_cell_and_msg(cell, exec_reply['content'])
nbclient.exceptions.CellExecutionError: An error occurred while executing the following cell:
------------------
import os
from sklearn.datasets import load_breast_cancer
from sklearn.tree import DecisionTreeClassifier
from sklearn.model_selection import train_test_split
from sklearn.metrics import confusion_matrix
from sklearn.tree import export_graphviz
from IPython.display import Image
from pydot import graph_from_dot_data
import pandas as pd
import numpy as np
cancer = load_breast_cancer()
X = pd.DataFrame(cancer.data, columns=cancer.feature_names)
print(X)
y = pd.Categorical.from_codes(cancer.target, cancer.target_names)
y = pd.get_dummies(y)
print(y)
X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=1)
tree_clf = DecisionTreeClassifier(max_depth=5)
tree_clf.fit(X_train, y_train)
export_graphviz(
tree_clf,
out_file="DataFiles/cancer.dot",
feature_names=cancer.feature_names,
class_names=cancer.target_names,
rounded=True,
filled=True
)
cmd = 'dot -Tpng DataFiles/cancer.dot -o DataFiles/cancer.png'
os.system(cmd)
------------------
---------------------------------------------------------------------------
ModuleNotFoundError Traceback (most recent call last)
<ipython-input-2-7c394b1e8b71> in <module>
 7 
 8 from IPython.display import Image
----> 9 from pydot import graph_from_dot_data
 10 import pandas as pd
 11 import numpy as np
ModuleNotFoundError: No module named 'pydot'
ModuleNotFoundError: No module named 'pydot'
@@ -0,0 +1,46 @@
Traceback (most recent call last):
File "/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/jupyter_cache/executors/utils.py", line 51, in single_nb_execution
executenb(
File "/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/nbclient/client.py", line 1087, in execute
return NotebookClient(nb=nb, resources=resources, km=km, **kwargs).execute()
File "/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/nbclient/util.py", line 74, in wrapped
return just_run(coro(*args, **kwargs))
File "/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/nbclient/util.py", line 53, in just_run
return loop.run_until_complete(coro)
File "/Users/mhjensen/opt/anaconda3/lib/python3.8/asyncio/base_events.py", line 616, in run_until_complete
return future.result()
File "/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/nbclient/client.py", line 540, in async_execute
await self.async_execute_cell(
File "/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/nbclient/client.py", line 832, in async_execute_cell
self._check_raise_for_error(cell, exec_reply)
File "/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/nbclient/client.py", line 740, in _check_raise_for_error
raise CellExecutionError.from_cell_and_msg(cell, exec_reply['content'])
nbclient.exceptions.CellExecutionError: An error occurred while executing the following cell:
------------------
heads_proba = 0.51
coin_tosses = (np.random.rand(10000, 10) < heads_proba).astype(np.int32)
cumulative_heads_ratio = np.cumsum(coin_tosses, axis=0) / np.arange(1, 10001).reshape(-1, 1)
plt.figure(figsize=(8,3.5))
plt.plot(cumulative_heads_ratio)
plt.plot([0, 10000], [0.51, 0.51], "k--", linewidth=2, label="51%")
plt.plot([0, 10000], [0.5, 0.5], "k-", label="50%")
plt.xlabel("Number of coin tosses")
plt.ylabel("Heads ratio")
plt.legend(loc="lower right")
plt.axis([0, 10000, 0.42, 0.58])
save_fig("votingsimple")
plt.show()
------------------
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
<ipython-input-1-eface79dac2c> in <module>
 1 heads_proba = 0.51
----> 2 coin_tosses = (np.random.rand(10000, 10) < heads_proba).astype(np.int32)
 3 cumulative_heads_ratio = np.cumsum(coin_tosses, axis=0) / np.arange(1, 10001).reshape(-1, 1)
 4 plt.figure(figsize=(8,3.5))
 5 plt.plot(cumulative_heads_ratio)
NameError: name 'np' is not defined
NameError: name 'np' is not defined
@@ -0,0 +1,31 @@
Traceback (most recent call last):
File "/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/jupyter_cache/executors/utils.py", line 51, in single_nb_execution
executenb(
File "/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/nbclient/client.py", line 1087, in execute
return NotebookClient(nb=nb, resources=resources, km=km, **kwargs).execute()
File "/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/nbclient/util.py", line 74, in wrapped
return just_run(coro(*args, **kwargs))
File "/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/nbclient/util.py", line 53, in just_run
return loop.run_until_complete(coro)
File "/Users/mhjensen/opt/anaconda3/lib/python3.8/asyncio/base_events.py", line 616, in run_until_complete
return future.result()
File "/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/nbclient/client.py", line 540, in async_execute
await self.async_execute_cell(
File "/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/nbclient/client.py", line 832, in async_execute_cell
self._check_raise_for_error(cell, exec_reply)
File "/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/nbclient/client.py", line 740, in _check_raise_for_error
raise CellExecutionError.from_cell_and_msg(cell, exec_reply['content'])
nbclient.exceptions.CellExecutionError: An error occurred while executing the following cell:
------------------
import numpy as np
x = np.log(np.array([4.0, 7.0, 8.0])
print(x)
------------------
 File "<ipython-input-6-f6d7a289d493>", line 3
 print(x)
 ^
SyntaxError: invalid syntax
SyntaxError: invalid syntax (<ipython-input-6-f6d7a289d493>, line 3)
@@ -0,0 +1,58 @@
Traceback (most recent call last):
File "/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/jupyter_cache/executors/utils.py", line 51, in single_nb_execution
executenb(
File "/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/nbclient/client.py", line 1087, in execute
return NotebookClient(nb=nb, resources=resources, km=km, **kwargs).execute()
File "/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/nbclient/util.py", line 74, in wrapped
return just_run(coro(*args, **kwargs))
File "/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/nbclient/util.py", line 53, in just_run
return loop.run_until_complete(coro)
File "/Users/mhjensen/opt/anaconda3/lib/python3.8/asyncio/base_events.py", line 616, in run_until_complete
return future.result()
File "/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/nbclient/client.py", line 540, in async_execute
await self.async_execute_cell(
File "/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/nbclient/client.py", line 832, in async_execute_cell
self._check_raise_for_error(cell, exec_reply)
File "/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/nbclient/client.py", line 740, in _check_raise_for_error
raise CellExecutionError.from_cell_and_msg(cell, exec_reply['content'])
nbclient.exceptions.CellExecutionError: An error occurred while executing the following cell:
------------------
# Blocking
@timeFunction
def blocking(self, blockSizeMax = 500):
blockSizeMin = 1
self.blockSizes = []
self.meanVec = []
self.varVec = []
for i in range(blockSizeMin, blockSizeMax):
if(len(self.data) % i != 0):
pass#continue
blockSize = i
meanTempVec = []
varTempVec = []
startPoint = 0
endPoint = blockSize
while endPoint <= len(self.data):
meanTempVec.append(np.average(self.data[startPoint:endPoint]))
startPoint = endPoint
endPoint += blockSize
mean, var = np.average(meanTempVec), np.var(meanTempVec)/len(meanTempVec)
self.meanVec.append(mean)
self.varVec.append(var)
self.blockSizes.append(blockSize)
self.blockingAvg = np.average(self.meanVec[-200:])
self.blockingVar = (np.average(self.varVec[-200:]))
self.blockingStd = np.sqrt(self.blockingVar)
------------------
 File "<ipython-input-6-2ff97f4bf03b>", line 2
 @timeFunction
 ^
IndentationError: unexpected indent
IndentationError: unexpected indent (<ipython-input-6-2ff97f4bf03b>, line 2)
File diff suppressed because one or more lines are too long
File diff suppressed because one or more lines are too long
@@ -1282,7 +1282,7 @@ C
K
$$
\frac{\partial \log{\vert\boldsymbol{A}\vert}}{\partial \boldsymbol{A}} = (\boldsymbol{A}^{-1})^T.
\frac{\partial\log{\vert\boldsymbol{A}\vert}}{\partial \boldsymbol{A}}=(\boldsymbol{A}^{-1})^T.
$$
The residuals $\boldsymbol{\epsilon}$ are in turn given by
@@ -2103,51 +2103,6 @@ $$
You can use the functionality included in scikit-learn. If you feel for it, you can use your own program and define functions which compute the above two functions.
Discuss the meaning of these results. Try also to vary the coefficient in front of the added stochastic noise term and discuss the quality of the fits.
!bsol
The code here is an example of where we define our own design matrix and fit parameters $\beta$.
import os
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split
def save_fig(fig_id):
plt.savefig(image_path(fig_id) + ".png", format='png')
def R2(y_data, y_model):
return 1 - np.sum((y_data - y_model) ** 2) / np.sum((y_data - np.mean(y_data)) ** 2)
def MSE(y_data,y_model):
n = np.size(y_model)
return np.sum((y_data-y_model)**2)/n
x = np.random.rand(100)
y = 2.0+5*x*x+0.1*np.random.randn(100)
# The design matrix now as function of a given polynomial
X = np.zeros((len(x),3))
X[:,0] = 1.0
X[:,1] = x
X[:,2] = x**2
# We split the data in test and training data
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)
# matrix inversion to find beta
beta = np.linalg.inv(X_train.T @ X_train) @ X_train.T @ y_train
print(beta)
# and then make the prediction
ytilde = X_train @ beta
print("Training R2")
print(R2(y_train,ytilde))
print("Training MSE")
print(MSE(y_train,ytilde))
ypredict = X_test @ beta
print("Test R2")
print(R2(y_test,ypredict))
print("Test MSE")
print(MSE(y_test,ypredict))
!esol
@@ -2229,7 +2184,5 @@ Perform an ordinary least squares and compute the means squared error and the $R
!bsubex
Add now a model which allows you to make polynomials up to degree $15$. Perform a standard OLS fitting of the training data and compute the MSE and $R2$ for the training and test data and plot both test and training data MSE and $R2$ as functions of the polynomial degree. Compare what you see with Figure 2.11 of Hastie et al. Comment your results. For which polynomial degree do you find an optimal MSE (smallest value)?
!bsol
Here you simply need to change the degree of the polynomial in the above code to $n=15$.
!esol
!esubex
Binary file not shown.

Before

Width:  |  Height:  |  Size: 6.1 KiB

After

Width:  |  Height:  |  Size: 6.2 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 16 KiB

Binary file not shown.

Before

Width:  |  Height:  |  Size: 9.3 KiB

After

Width:  |  Height:  |  Size: 9.4 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 21 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 25 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 95 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 33 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 10 KiB

Binary file not shown.

Before

Width:  |  Height:  |  Size: 14 KiB

After

Width:  |  Height:  |  Size: 13 KiB

Binary file not shown.

Before

Width:  |  Height:  |  Size: 9.8 KiB

After

Width:  |  Height:  |  Size: 9.9 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 16 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 19 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 26 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 25 KiB

+2 -64
View File
@@ -1963,7 +1963,7 @@
"metadata": {},
"source": [
"$$\n",
"\\frac{\\partial \\log{\\vert\\boldsymbol{A}\\vert}}{\\partial \\boldsymbol{A}} = (\\boldsymbol{A}^{-1})^T.\n",
"\\frac{\\partial\\log{\\vert\\boldsymbol{A}\\vert}}{\\partial \\boldsymbol{A}}=(\\boldsymbol{A}^{-1})^T.\n",
"$$"
]
},
@@ -3359,66 +3359,6 @@
"You can use the functionality included in scikit-learn. If you feel for it, you can use your own program and define functions which compute the above two functions. \n",
"Discuss the meaning of these results. Try also to vary the coefficient in front of the added stochastic noise term and discuss the quality of the fits.\n",
"\n",
"!bsol\n",
"The code here is an example of where we define our own design matrix and fit parameters $\\beta$."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false,
"editable": true
},
"outputs": [],
"source": [
"import os\n",
"import numpy as np\n",
"import pandas as pd\n",
"import matplotlib.pyplot as plt\n",
"from sklearn.model_selection import train_test_split\n",
"\n",
"def save_fig(fig_id):\n",
" plt.savefig(image_path(fig_id) + \".png\", format='png')\n",
"\n",
"def R2(y_data, y_model):\n",
" return 1 - np.sum((y_data - y_model) ** 2) / np.sum((y_data - np.mean(y_data)) ** 2)\n",
"def MSE(y_data,y_model):\n",
" n = np.size(y_model)\n",
" return np.sum((y_data-y_model)**2)/n\n",
"\n",
"x = np.random.rand(100)\n",
"y = 2.0+5*x*x+0.1*np.random.randn(100)\n",
"\n",
"\n",
"# The design matrix now as function of a given polynomial\n",
"X = np.zeros((len(x),3))\n",
"X[:,0] = 1.0\n",
"X[:,1] = x\n",
"X[:,2] = x**2\n",
"# We split the data in test and training data\n",
"X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)\n",
"# matrix inversion to find beta\n",
"beta = np.linalg.inv(X_train.T @ X_train) @ X_train.T @ y_train\n",
"print(beta)\n",
"# and then make the prediction\n",
"ytilde = X_train @ beta\n",
"print(\"Training R2\")\n",
"print(R2(y_train,ytilde))\n",
"print(\"Training MSE\")\n",
"print(MSE(y_train,ytilde))\n",
"ypredict = X_test @ beta\n",
"print(\"Test R2\")\n",
"print(R2(y_test,ypredict))\n",
"print(\"Test MSE\")\n",
"print(MSE(y_test,ypredict))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"!esol\n",
"\n",
"\n",
"\n",
@@ -3545,9 +3485,7 @@
"\n",
"!bsubex\n",
"Add now a model which allows you to make polynomials up to degree $15$. Perform a standard OLS fitting of the training data and compute the MSE and $R2$ for the training and test data and plot both test and training data MSE and $R2$ as functions of the polynomial degree. Compare what you see with Figure 2.11 of Hastie et al. Comment your results. For which polynomial degree do you find an optimal MSE (smallest value)?\n",
"!bsol\n",
"Here you simply need to change the degree of the polynomial in the above code to $n=15$.\n",
"!esol\n",
"\n",
"!esubex"
]
}