• Write a program to estimate the Euler constant by Monte Carlo simulation.

(1)

Exercise: e ∼ 2.7183

• Write a program to estimate the Euler constant by Monte Carlo simulation.

• It can be done as follows.

• Let N be the number of iterations.

• For each iteration, find the minimal number n so that P

n

i =1

r

_i

> 1 where r

_i

is the random variable following the standard uniform distribution (you can simply use rand).

• Then e is the average of n.

(2)

Special Issue: Sort

1 >> stocks = {"GOOG", 15;

2 "TSMC", 12;

3 "AAPL", 18};

4 >> [~, idx] = sort([stocks{:, 2}], "descend")

5

6 idx =

7

8 3 1 2

9

10 >> stocks = stocks(idx, :)

11

12 stocks =

13

14 "AAPL" [18]

15 "GOOG" [15]

16 "TSMC" [12]

Zheng-Liang Lu 101

(3)

Programming Exercise: Sorting Algorithm

¹

• Let A be any array.

• Write a program which outputs the sorted array of A (in ascending order).

• For example, A = [5, 4, 1, 2, 3].

• Then the sorted array is [1, 2, 3, 4, 5].

(4)

Special Issue: Random Permutation

• Use randperm to generate an index array with a random order.

1 >> A = ["Matlab", "Python", "Java", "C++"];

2 >> idx = randperm(length(A))

3

4 idx =

5

6 3 1 2 4

7

8 >> A(idx)

9

10 ans =

11

12 1x4 string array

13

14 "Java" "Matlab" "Python" "C++"

Zheng-Liang Lu 103

(5)

“Exploring the unknown requires tolerating uncertainty.”

– Brian Greene

“I can live with doubt, and uncertainty, and not knowing.

I think it is much more interesting to live not knowing than have answers which might be wrong.”

– Richard Feynman

(6)

Speedup: Vectorization (Revisited)

²

• Vector in, vector out.

1 >> x = randi(100, 1, 5)

2

3 x =

4

5 88 30 90 73 82

6

7 >> dx = diff(x)

8 9 dx =

10

11 -58 60 -17 9

2More aboutvectorization.

Zheng-Liang Lu 105

(7)

Advantages from Vectorization

• Appearance: vectorized mathematical code appears more like the mathematical expressions found in textbooks, making the code easier to understand.

• Less error prone: without loops, vectorized code is often shorter.

• Fewer lines of code mean fewer opportunities to introduce programming errors.

• Performance: vectorized code often runs much faster than the

corresponding code containing loops.

(8)

Performance Analysis: Profiling

• Use a timer to measure your performance.

³

• In newer version, press the button Run and Time.

• Identify which functions are consuming the most time.

• Know why you are calling them and then look for alternatives to improve the overall performance.

3Note that the results may differ depending on the difference of run-time environments, so make sure that you benchmark the algorithms on thesame conditions.

Zheng-Liang Lu 107

(9)

tic & toc

• The command tic makes a stopwatch timer start.

• The command toc returns the elapsed time from the stopwatch timer started by tic.

1 >> tic

2 >> toc

3 Elapsed time is 0.786635 seconds.

4 >> toc

6 >> toc

(10)

Selected Performance Suggestions

⁴

• Preallocate arrays.

• Instead of continuously resizing arrays, consider preallocating the maximum amount of space required for an array.

• Vectorize your code.

• Create new variables if data type changes.

• Use functions instead of scripts.

• Avoid overloading Matlab built-in functions.

4SeeTechniques for Improving Performance.

Zheng-Liang Lu 109

(11)

Programming Exercise: A Benchmark

• Let N = 1e1, 1e2, 1e3, 1e4, 1e5.

• Write a program which produces a benchmark for the following three cases:

• Generate an array of 1 : N by dynamically resizing the array.

• Generate an array of 1 : N by allocating an array of size N and filling up sequentially.

• Generate an array of 1 : N by vectorization.

(12)

Analysis of Algorithms (Optional)

• For one problem, there exist various algorithms (solutions).

• We then compare these algorithms for various considerations and choose the most appropriate one.

• In general, we want efficient algorithms.

• Except for real-time performance analysis, could we predict before the program is completed?

• Definitely yes.

Zheng-Liang Lu 111

(13)

Growth Rate

• Now we use f (n) to denote the growth rate of time cost as a function of n.

• In general, n refers to the data size.

• For simplicity, assume that every instruction (e.g. + − ×÷) takes 1 unit of computation time.

• Find f (n) for the following problem.

• Sum(n): ?

• Triangle(n): ?

(14)

O-notation

⁵

• In math, O-notation describes the limiting behavior of a function, usually in terms of simple functions.

• We say that

f (n) ∈ O(g (n)) as n → ∞ if and only if ∃c > 0, n

₀

> 0 such that

| f (n) | ≤ c| g (n) | ∀n ≥ n

₀

.

• So O(g (n)) is a collection featured by a simple function g (n).

• We use f (n) ∈ O(g (n)) to denote that f (n) is one instance of O(g (n)).

5See https://en.wikipedia.org/wiki/Big_O_notation.

Zheng-Liang Lu 113

(15)

• Big-O is used for the asymptotic upper bound of time

complexity of algorithm.

(16)

• For example, 8n

²

− 3n + 4 ∈ O(n

²

).

• For large n, you could ignore the last two terms. (Why?)

• It is easy to find a constant c > 0 so that cn²> 8n², say c = 9.

• Hence the statement is proved.

• Also, 8n

²

− 3n + 4 ∈ O(n

³

) but we seldom say this. (Why?)

• However, 8n

²

− 3n + 4 / ∈ O(n). (Why?)

• What is this analysis related to the algorithm?

• Any insight?

Zheng-Liang Lu 115

(17)

Common Simple Functions

⁶

(18)

Remarks

• We often make a trade-off between time and space.

• Unlike time, we can reuse memory.

• Users are sensitive to time.

• Playing game well is hard.

⁷

• Solve the problem P ?= NP, which is one of Millennium Prize Problems.

⁸

7See https://en.wikipedia.org/wiki/Game_complexity.

8See https://en.wikipedia.org/wiki/P_versus_NP_problem.

Zheng-Liang Lu 117

(19)

“All roads lead to Rome.”

– Anonymous

“但如你根本並無招式，敵人如何來破你的招式？”

– 風清揚。笑傲江湖。第十回。傳劍

(20)

1 >> Lecture 3

2 >>

3 >> -- Graphics

4 >>

Zheng-Liang Lu 119

(21)

Introduction

• Engineers use graphic techniques to make the information easier to understand.

• With graphs, it is easy to identify trends, pick out highs and lows, and isolate data points that may be measurement or calculation errors.

• Graphs can also be used as a quick check to determine if a computer solution is yielding expected results.

• A set of ordered pairs is used to identify points on a 2D graph.

(22)

2D Line Plot

• plot(x , y ) creates a 2D line plot for all (x , y ) pairs in order.

• You may use more parameters for the plot as follows:

Zheng-Liang Lu 121

(23)

Example

1 clear; clc; close all;

2

3 x = linspace(0, 2 * pi, 20);

4 y = sin(x);

5

6 figure; plot(x, y, "r-o");

7 grid on;

• Call figure to create a figure.

• Use close to close all figures or specific one.

• Use grid to add the gray grid as the background.

(24)

0 1 2 3 4 5 6 7 -1

-0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1

Zheng-Liang Lu 123

(25)

Example: Multiple Curves

2

3 x = linspace(0, 2 * pi, 50);

4 figure; hold on; grid on;

5 plot(x, sin(x), '*');

6 plot(x, cos(x), 'o');

7 plot(x, sin(x) + cos(x), '+');

• Use hold to put multiple curves in the same figure.

(26)

0 1 2 3 4 5 6 7 -1.5

-1 -0.5 0 0.5 1 1.5

Zheng-Liang Lu 125

(27)

Selected Annotations

• Use title to add a title to the plot.

• Use xlabel to add a label to the x axis of the plot.

• Use ylabel to add a label to the y axis of the plot.

• Use legend to add legends for lines.

• More annotations can be created by annotation.

⁹

• Note that you can always generate the codes associated with

the plot you modified.

(28)

Example

2

3 x = linspace(0, 2 * pi, 30);

4 y = sin(x); z = cos(x);

5

7 plot(x, y, "o--");

8 plot(x, z, "x:");

9 legend("Input", "Output", "location", "best");

10

11 xlabel("Time (s)"); ylabel("Amplitude (unit)");

12 title("Title");

13 annotation("textarrow", [.3, .6], [.7, .4] , ...

14 "String", "ABC");

Zheng-Liang Lu 127

(29)

-0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1

Amplitude (unit)

Title

Input Output ABC

(30)

Graphics Objects

• You can use plot tool (in the figures) to change the properties.

• Graphics objects are the components for data visualization.

• Each object can be assigned to a unique identifier, called a graphics handle.

• Via graphics handles, you can manipulate their properties

¹⁰

by the following instructions:

• set: set properties.

• get: query properties.

10See http:

//www.mathworks.com/help/matlab/graphics-object-properties.html.

Zheng-Liang Lu 129

(31)

Example

2

3 x = linspace(-1, 1, 100);

4 h = plot(x, sin(1 ./ x));

5 grid on;

6 set(h, "Marker", "o");

7 set(h, "MarkerSize", 5);

8 set(h, "LineWidth", 1.5);

9 set(h, "LineStyle", ":");

10 set(h, "MarkerEdgeColor", "r");

11 set(h, "MarkerFaceColor", "g");

(32)

-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 -1

-0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1

Zheng-Liang Lu 131

(33)

Graphics Object Identification

¹¹

• gcf: get current figure

• gca: get current axis

• gco: get current object

(34)

Output Figures

• You can save one figure as a specific image format.

• For example, bmp, jpeg, and eps.

• Use the hot key ctrl + s .

Zheng-Liang Lu 133

(35)

• You can also use print to save the figures.

¹²

2

3 x = linspace(0, 2 * pi, 20);

4 y = sin(x);

5

6 figure; plot(x, y, "r-o"); grid on;

7 print(gcf, "-djpeg", "sin.jpg", "-r300");

• Use saveas to save figure in a specific file format.

¹³

• Use savefig to save figure and contents to fig-file.

¹⁴

12

(36)

Exercise: TWSE:IND

2

3 [~, ~, raw] = xlsread("y9999.xlsx");

4 prices = [raw{4 : end, 2}];

5 volumes = [raw{4 : end, 3}];

6 dates = datetime(raw(4 : end, 1), ...

7 "format", "yyyy/MM/dd");

8

9 fig1 = figure;

10 plot(dates, prices); grid on;

11 ylabel("TWSE:IND");

12 annotation(fig1, "arrow", [0.4 0.88], [0.28 0.65]);

• Use datetime to convert a date string to a datetime object.

• Note that you need to specify a date format, say

"yyyy/MM/dd"

.

Zheng-Liang Lu 135

(37)

5000 6000 7000 8000 9000 10000 11000 12000 13000

TWSE:IND

(38)

Bar Plot

¹⁵

• Use bar to draws a bar chart, for example,

2

3 x = randi(100, 8, 3);

4 bar(x); grid on;

• Try barh.

15See http://www.mathworks.com/help/matlab/ref/bar.html and http://www.mathworks.com/help/matlab/creating_plots/

overlay-bar-graphs.html.

Zheng-Liang Lu 137

(39)

10 20 30 40 50 60 70 80 90 100

(40)

Exercise: Traded Volumes of TWSE:IND

2

8

9 figure;

10 bar(dates, volumes); grid on;

11 ylabel("Traded volume");

Zheng-Liang Lu 139

(41)

2 4 6 8 10 12

Traded volume

10⁶

(42)

Dual y-Axes Plot

• Use yyaxis to specify the left/right y axis, for example,

2

3 x = linspace(0, 30, 300);

4 y1 = 10 * exp(-0.05 * x) .* sin(x);

5 y2 = 0.8 * exp(-0.5 * x) .* sin(10 * x);

6

7 figure;

8 yyaxis left;

9 plot(x, y1); grid on;

10 yyaxis right;

11 plot(x, y2);

• Use plotyy in old version.

Zheng-Liang Lu 141

(43)

-6 -4 -2 0 2 4 6 8 10

-0.4 -0.2 0 0.2 0.4 0.6 0.8

(44)

Exercise: Index feat. Volume in One Figure

2

8

9 yyaxis left; plot(dates, prices);

10 ylabel("TWSE:IND"); grid on;

11 yyaxis right; bar(dates, volumes);

12 ylabel("Traded volume"); grid on;

Zheng-Liang Lu 143

(45)

5000 6000 7000 8000 9000 10000 11000 12000 13000

TWSE:IND

2 4 6 8 10 12

Traded volume

10⁶

(46)

Histogram Plot

¹⁷

• Histograms group the numeric data into bins.

• Use histogram to create histogram plots.

¹⁶

2

3 data = randn(1, 1e3) .ˆ 2;

4 figure;

5 histogram(data, ...

6 "BinMethod", "integers", ...

7 "Normalization", "probability");

8 grid on;

16If your version is before 2014, use hist.

17More details could be found in https://www.mathworks.com/help/

matlab/ref/matlab.graphics.chart.primitive.histogram.html.

Zheng-Liang Lu 145

(47)

0.1 0.2 0.3 0.4 0.5 0.6

(48)

Exercise: Distribution of Return Rates of TWSE:IND

2

8 return rates = diff(prices) ./ prices(1 : end - 1);

9

10 figure;

11 histogram(return rates * 100, ...

12 "binmethod", "integer", ...

13 "normalization", "probability");

14 xlabel("Return rate (%)");

15 ylabel("Probability"); grid on;

Zheng-Liang Lu 147

(49)

-8 -6 -4 -2 0 2 4 6 8 0

0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5

Probability

(50)

Grid Plot: subplot

¹⁸

2008 2010 2012 2014 2016 2018 2020

6000 8000 10000 12000

TWSE:IND

20060 2008 2010 2012 2014 2016 2018 2020

5 10

Traded volume

10⁶

18See https://www.mathworks.com/help/matlab/ref/subplot.html.

Zheng-Liang Lu 149

(51)

2

8

9 figure;

10 subplot(2, 1, 1); plot(dates, prices); grid on;

11 ylabel("TWSE:IND");

12 subplot(2, 1, 2); bar(dates, volumes, "r"); grid on;

13 ylabel("Traded volume");

• Use subplot(m, n, p) to divide the current figure into an

m-by-n grid and use p to specify the certain subplot.

(52)

Digression: Table

¹⁹

• Use table to create a table for column-oriented or tabular data that is often stored as columns in a spreadsheet.

• Use detectImportOptions to create import options based on the contents of a file (if readtable cannot read files correctly).

• Use stackedplot to draw a stacked plot of several variables with common x-axis.

19See https://www.mathworks.com/help/matlab/tables.html.

Zheng-Liang Lu 151

(53)

2

3 filename = "2330.xlsx";

4 s2330 = readtable(filename, ...

5 detectImportOptions(filename));

6 % Delete the first two rows.

7 s2330(1 : 2, :) = [];

8 % Assign the header name for each column.

9 s2330.Properties.VariableNames = ["Date", "Open", ...

"High", "Low", "Close", "Volume"];

10 % Convert date strings to datetime objects.

11 s2330.Date = datetime(s2330.Date, ...

12 "format", "yyyy-MM-dd");

13 % Use stackedplot to draw an interactive plot!

14 stackedplot(s2330, {"Close", "Volume"}, ...

15 "xvar", "Date"); grid on;

(54)

200 250 300

Close

Jul 2017 Jan 2018 Jul 2018 Jan 2019 Jul 2019 Date

2 4 6 8 10 12 14

Volume 10⁴

Zheng-Liang Lu 153

(55)

Selected Table Functions

• File I/O: readtable, writetable.

• Summary information: head, tail, summary, stackedplot.

• Sort, rearrange, and customize: sortrows, unique, addvars, removevars, rows2vars, stack, unstack, inner2outer.

• Join and set operations: join, innerjoin, outerjoin, union, intersect, ismember, setdiff, setxor.

• Apply functions to table contents: varfun, rowfun,

findgroups, splitapply, groupsummary

(56)

Exercise: Merging Two Tables

2

3 gspc = readtable("ˆGSPC.csv");

4 twii = readtable("ˆTWII.csv", ...

5 detectImportOptions("ˆTWII.csv"));

6 twii.Date = datetime(twii.Date, ...

7 "format", gspc.Date.Format);

8 % Merge two time series by union of dates.

9 merged table = outerjoin(twii, gspc, ...

10 "Keys", "Date", ...

11 "MergeKeys", 1);

12 stackedplot(merged table, ...

13 {"Close twii", "Close gspc"}, ...

14 "xvariable", "Date"); grid on;

Zheng-Liang Lu 155

(57)

8000 9000 10000 11000

Close_twii

2000 2200 2400 2600 2800

Close_gspc

(58)

Candle Chart with Timetable

²⁰

1 % Ingore the part identical to the previous ...

example of 2330.

2

3 s2330 = table2timetable(s2330, "RowTimes", "Date");

4 candle(s2330(end - 30 : end, :)); % last 30 days

5 ylabel("Daily price (TWD)");

• Use timetable to convert the table (with variable names:

"Open"

,

"High"

,

"Low"

,

"Close"

) to a timetable by specifying the RowTimes.

• Try priceandvol.

20See https://www.mathworks.com/help/finance/candle.html and https://www.mathworks.com/help/matlab/timetables.html with https://www.mathworks.com/help/finance/examples/

using-timetables-in-finance.html.

Zheng-Liang Lu 157

(59)

300 305 310 315 320 325 330 335 340 345

Daily price (TWD)

(60)

Smart Plot: fplot

• Use fplot to make a line plot over a specific range with adaptive steps.

• You may assign a function in a string form to fplot.

²¹

2

3 fplot("sin(1 / x)", [-0.05, 0.05]); grid on;

21Warning: fplot will not accept character vector or string inputs in a future release.

Zheng-Liang Lu 159

(61)

-0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1

(62)

errorbar

²²

2

3 x = 0 : 10 : 100;

4 y = [20 30 45 40 60 65 80 75 95 90];

5 yneg = [1 3 5 3 5 3 6 4 3 3];

6 ypos = [2 5 3 5 2 5 2 2 5 5];

7 xneg = [1 3 5 3 5 3 6 4 3 3];

8 xpos = [2 5 3 5 2 5 2 2 5 5];

9 errorbar(x, y, yneg, ypos, xneg, xpos, "o");

10 grid on;

22See https://www.mathworks.com/help/matlab/ref/errorbar.html.

Zheng-Liang Lu 161

(63)

20 30 40 50 60 70 80 90 100

(64)

Pie Chart

26%

19%

21%

13%

21%

Zheng-Liang Lu 163

(65)

2

3 X = rand(1, 5);

4 labels = {"A", "B", "C", "D", "E"};

5 explode = [0, 1, 0, 1, 0];

6 pie(X, explode, labels);

• Use pie to create a pie chart.

²³

• Note that the explode vector is used to offset slices for the

nonzero elements.

(66)

Word Cloud

market_values

Zheng-Liang Lu 165

(67)

2

3 [~, ~, raw] = xlsread("twse mktValue.xlsx");

4

5 stock ticks = string(raw(4 : end, 1));

6 idx = strcmp(raw(:, 3), "-"); % Find all "-"s.

7 raw(idx, 3) = {0}; % Replace them by 0.

8 market values = [raw{4 : end, 3}]';

9

10 tbl = table(stock ticks, market values);

11 figure;

12 wordcloud(tbl, "stock ticks", "market values");

• Use strcmp to compare strings and return a boolean vector.

• Use wordcloud to create a word cloud chart from text data.

²⁴

(68)

Contours

²⁵

• Use meshgrid to partition the specified range of x and y .

• Note that the return values are in form of matrices. (Why?)

2

3 x = linspace(-2 * pi, 2 * pi);

4 y = linspace(0, 4 * pi);

5 [X, Y] = meshgrid(x, y);

6 Z = sin(X) + cos(Y); % Using vectorization.

7 figure; contour(X, Y, Z, "showtext", "on");

25See https://www.mathworks.com/help/matlab/ref/contour.html.

You may try contourf.

Zheng-Liang Lu 167

(69)

-1.5 -1.5

-1

-1 -1

-0.5 -0.5 -0.5

-0.5

-0.5 -0.5

-0.5

0 0

0

0 0

0

0 0

0.5 0.5

0.5

0.5 0.5

1 1

1

1.5 1.5

2 4 6 8 10 12

(70)

Quiver (Velocity) Plot

• Use quiver(x , y , u, v ) to plot a vector (u, v ) at the coordinate (x , y ).

• Use peaks with a positive number as sample size to generate a set of 3d points.

²⁶

2

3 [x, y, z] = peaks(20);

4 [u, v] = gradient(z);

6 contour(x, y, z, 10);

7 quiver(x, y, u, v);

26See https://www.mathworks.com/help/matlab/ref/peaks.html.

Zheng-Liang Lu 169

(71)

-2 -1 0 1 2 3

(72)

Mesh Plot

• Use mesh to draw a wireframe mesh.

2

3 x = linspace(-3, 3, 50);

4 y = x + pi / 2;

5 [X, Y] = meshgrid(x, y);

6 Z = cos(X) .* sin(Y);

7 figure; mesh(X, Y, Z); grid on;

8 xlabel("x"); ylabel("y"); zlabel("z");

Zheng-Liang Lu 171

(73)

-1 6 -0.5

4 4

z 0

2 0.5

2 1

0 0

(74)

Surface Plot

• Use surf to draw a colored surface.

• Try meshz, meshc, surfc, and waterfall.

2

3 x = linspace(-2, 2, 25);

4 y = linspace(-2, 2, 25);

5 [X, Y] = meshgrid(x, y); % form all x-y pairs

6 Z = X .* exp(-X .ˆ 2 - Y .ˆ 2);

7 surf(X, Y, Z);

8 xlabel("x"); ylabel("y"); zlabel("z");

Zheng-Liang Lu 173

(75)

-0.5 2

1 2

z 0

0 1 0.5

-1 0

(76)

Exercise

• Write a program to draw a surface plot for sinc(R) =

^sin(R)_R

.

• Note that there exists a singularity at R = 0, which should be removed by replacing a zero with eps = 2.2204 × 10

⁻¹⁶

.

2

3 [X, Y] = meshgrid(linspace(-10, 10, 51));

4 R = sqrt(X .ˆ 2 + Y .ˆ 2);

5 R(R == 0) = eps; % Avoid the singularity.

6 Z = sin(R) ./ R;

7 surf(X, Y, Z);

Zheng-Liang Lu 175

(77)

-0.5 10 0

5 10

0.5

0 5 1

-5 0

(78)

3D Line Plot

• Use plot3 to draw a 3d curve.

2

3 t = linspace(0, 10 * pi, 100);

4 x = t .* sin(t); y = t .* cos(t);

5

6 figure;

7 plot3(x, y, t); hold on;

8 plot3(x, y, -t, "r");

9 axis equal; grid on;

Zheng-Liang Lu 177

(79)

-30 -20 -10

20 0

20 10

20

0 0

30

(80)

Misc

²⁹

• Use the button Rotate 3D to change the view angle.

• Use view to set the view angle.

²⁷

• Try colorbar and colormap.

²⁸

2

3 peaks;

4 view([117, 58]); % View angle in degree.

5 colorbar; % Appends a colorbar to the current axes.

6 colormap summer; % Change the colormap.

27az: azimuth (horizontal) rotation; el: vertical elevation.

28See https://www.mathworks.com/help/matlab/ref/colormap.html.

29See

https://www.mathworks.com/products/matlab/plot-gallery.html.

Zheng-Liang Lu 179

(81)

-3 -2 -1 -3 0

-2 -5

-1 1

0 2

Peaks

0 5

-6 -4 -2 0 2 4 6 8