851 0585 04 L Modeling and Simulating Social

851 -0585 -04 L – Modeling and Simulating Social Systems with MATLAB Lecture 2 – Statistics and Plotting in MATLAB Karsten Donnay and Stefano Balietti Chair of Sociology, in particular of Modeling and Simulation © ETH Zürich | 2011 -10 -03

Goals of Lecture 2: students will 1. Consolidate knowledge of basic MATLAB functionalities acquired during lecture 1 through brief repetition of the main points and revision of the exercises. 2. Get familiar with some basic statistical functionality of MATLAB. 3. Understand the basic MATLAB 2 D plotting functionality and apply it appropriately; the example of logarithmic plots is discussed in detail. 4. Master positioning multiple graphs in the screen with subplot() and figure() and custom. m files. 2011 -10 -03 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 2

Repetition § Creating a row vector: >> v=[1 2 3 4 5] v = 1 2 3 4 5 2011 -10 -03 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 3

Repetition § Creating a row vector, 2 nd method: >> v=1: 5 v = 1 2 3 4 5 2011 -10 -03 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 4

Repetition § Creating a row vector, 3 rd method: linspace(start. Val, end. Val, n) >> v=linspace(0. 5, 1. 5, 5) v = 0. 50 0. 75 1. 00 1. 25 1. 50 2011 -10 -03 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 5

Repetition § Creating a column vector: >> v=[1; 2; 3; 4; 5] v = 1 2 3 4 5 2011 -10 -03 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 6

Repetition § Creating a column vector, 2 nd method: >> v=[1 2 3 4 5]’ v = 1 2 3 4 5 2011 -10 -03 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 7

Repetition § Creating a matrix: >> A=[1 2 3 4; 5 6 7 8] A = 1 2 3 4 5 6 7 8 2011 -10 -03 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 8

Repetition § Creating a matrix, 2 nd method: >> A=[1: 4; 5: 8] A = 1 2 3 4 5 6 7 8 2011 -10 -03 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 9

Repetition § Accessing an element in a vector: >> v = 1: 10; >> v(2) ans = 2 2011 -10 -03 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 10

Repetition § Accessing an element in a vector: >> v = 1: 10; >> v(2) ans = Vector indexes are 1 -based index = 1. . size ( v , 2 ) 2 1: column vector 2: row vector 2011 -10 -03 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 11

Repetition § Accessing an element in a matrix: A(row. Index, column. Index) >> A(2, 1) ans = 5 2011 -10 -03 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 12

Repetition - operators § Scalar operators: Basic: +, -, *, / Exponentialisation: ^ Square root: sqrt() 2011 -10 -03 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 13

Repetition - operators § Matrix operators: Basic: +, -, * § Element-wise operators: Multiplication: . * Division: . / Exponentialisation: . ^ § Solving Ax=b: x = Ab 2011 -10 -03 Attention! / (slash) and (back slash) produce different results K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 14

Repetition – for loop § Computation can be automized with for loops: >> y=0; for x=1: 4 y = y + x^2 + x; end >> y y = 40 2011 -10 -03 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 15

Repetition – if case § Conditional computation can be made with if: >> val=-4; if (val>0 ) abs. Val = val; else abs. Val = -val; end 2011 -10 -03 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 16

Lesson 1: Exercise 1 § Compute: a) b) c) 2011 -10 -03 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 17

Lesson 1: Exercise 1 – solution § Compute: a) >> (18+107)/(5*25) ans = 1 2011 -10 -03 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 18

Lesson 1: Exercise 1 – solution § Compute: b) >> s=sum(1: 100) or >> s=sum(1: 1: 100) default value or >> s=sum(linspace(1, 100)) s = default value 5050 2011 -10 -03 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 19

Lesson 1: Exercise 1 – solution § Compute: c) >> s=0; >> for i=5: 10 >> s=s+i^2 -i; >> end >> s s = 310 2011 -10 -03 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 20

Lesson 1: Exercise 2 § Solve for x: 2011 -10 -03 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 21

Lesson 1: Exercise 2 – solution >> A=[2 -3 -1 4; 2 3 -3 2; 2 -1 -1 -1; 2 -1 2 5]; >> b=[1; 2; 3; 4]; >> x=Ab x = 1. 9755 0. 3627 0. 8431 -0. 2549 Ax=b >> A*x ans = 1. 0000 2. 0000 3. 0000 4. 0000 2011 -10 -03 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 22

Lesson 1: Exercise 3 § Fibonacci sequence: Write a function which computes the Fibonacci sequence until a given number n and return the result in a vector. § The Fibonacci sequence F(n) is given by : 2011 -10 -03 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 23

Lesson 1: Exercise 3 – iterative solution fibonacci. m: function [v] = Fibonacci(n) v(1) = 0; if ( n>=1 ) v(2) = 1; end for i=3: n+1 v(i) = v(i-1) + v(i-2); end >> Fibonacci(7) ans = 0 1 1 2 3 5 8 13 2011 -10 -03 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 24

Lesson 1: Exercise 3 – recursive solution fibo_rec. m: function f = fibo_rec(n) if n == 0 f(1) = 0; elseif n == 1 f(2) = 1; elseif n > 1 f = fibo_rec(n - 1); f(n + 1) = f(n) + f(n - 1); end >> fibo_rec(7) ans = 0 1 1 2 3 5 8 13 2011 -10 -03 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 25

Statistics functions § Statistics in MATLAB can be performed with the following commands: Vector-based Mean value: mean(x) Median value: median(x) Min/max values: min(x), max(x) Standard deviation: std(x) Variance: var(x) Matrix-based Covariance: cov(x) Correlation coefficient: corrcoef(x) 2011 -10 -03 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 26

Where can I get help? >>help functionname § >>doc functionname § >>helpwin § Click on “More Help” from contextual pane after opening first parenthesis, e. g. : plot(… § Click on the Fx symbol before command prompt. 2011 -10 -03 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 27

2 D Plotting functions § 2 D plotting of curves (y vs x) in MATLAB Linear scale: plot(x, y) Double-logarithmic scale: loglog(x, y) Semi-logarithmic scale: semilogx(x, y) semilogy(x, y) Plot histogram of x: hist(x) 2011 -10 -03 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 28

Logarithmic Scale § It uses the logarithm of a physical quantity instead of the quantity itself. § E. g. , increment on the axis are powers of 10 (1, 100, 1000, …) instead of 1, 2, 3, 4. § Examples: 2011 -10 -03 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 29

Logarithmic Scale § It uses the logarithm of a physical quantity instead of the quantity itself. § E. g. , increment on the axis are powers of 10 (1, 100, 1000, …) instead of 1, 2, 3, 4. § Examples: Same distance between 1 -2 -3 and 10 -20 -30 2011 -10 -03 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 30

Logarithmic Rescaling: example ? 2011 -10 -03 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 31

Logarithmic Rescaling: example 2011 -10 -03 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 32

Logarithmic Scale: examples § Decibel § Richter scale § Entropy in thermodynamics § Information in information theory. § PH for acidity and alkalinity; § Stellar magnitude scale for brightness of stars § bit [log 2] and the byte 8[log 2] 2011 -10 -03 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 33

Logarithmic Scale: why and when to use it § If the data cover a large range of values the logarithm reduces this to a more manageable range. § Some of our senses operate in a logarithmic fashion (Weber–Fechner law) § Some types of equations get simplified: § Exponential: Y = ea. X -> semilog scale -> Y = -a. X § Power Law: Y = Xb -> log-log scale -> log. Y = b log. X 2011 -10 -03 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 34

Logarithmic Scale 2011 -10 -03 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 35

Back to plotting… details of plot § An additional parameter can be provided to plot() to define how the curve will look like: plot(x, y, ‘key’) Where key is a string which can contain: Color codes: ‘r’, ‘g’, ‘b’, ‘k’, ‘y’, … Line codes: ‘-’, ‘--’, ‘. -’ (solid, dashed, etc. ) Marker codes: ‘*’, ‘s’, ’x’ Examples: plot(x, y, ‘r--’) plot(x, y, ‘g*’) * * 2011 -10 -03 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 36

Plotting tips § To make the plots look nicer, the following commands can be used: Set label on x axis: xlabel(‘text’) Set label on y axis: ylabel(‘text’) Set title: title(‘text’) 2011 -10 -03 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 37

Plotting tips § Two additional useful commands: hold on|off § grid on|off § >> x=[-5: 0. 1: 5]; >> y 1=exp(-x. ^2); >> y 2=2*exp(-x. ^2); >> y 3=exp(-(x. ^2)/3); 2011 -10 -03 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 38

Plotting tips >> plot(x, y 1); >> hold on >> plot(x, y 2, ’r’); >> plot(x, y 3, ’g’); 2011 -10 -03 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 39

Plotting tips >> plot(x, y 1); >> hold on >> plot(x, y 2, ’r’); >> plot(x, y 3, ’g’); >> grid on 2011 -10 -03 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 40

Subplot() and figure() § subplot(m, n, p) breaks the figure window into an m-by -n matrix and selects the pth axes object for the current plot. § The axes are counted along the top row of the figure window, then the second row, etc. § h = figure(. . . ) adds a new figure window to the screen. It also returns the handle to it, so it can be further processed. § E. g. : saveas(h, ’myimage. jpg’); 2011 -10 -03 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 41

Subplot() by the way… We have already seen an example of a chart created using subplot on page 32. Could you reproduce it? 2011 -10 -03 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 42

Exercises: working with datasets § Two datasets for statistical plotting can be found on the course web page (with NETHZ login) http: //www. soms. ethz. ch/matlab you will find the files: countries. m cities. m 2011 -10 -03 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 43

Exercises: Datasets § Download the files countries. m and cities. m and save them in the working directory of MATLAB. 2011 -10 -03 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 44

Exercises: Datasets – countries § This dataset countries. m contains a matrix A with the following specification: § Rows: Different countries § Column 1: Population § Column 2: Annual growth (%) § Column 3: Percentage of youth § Column 4: Life expectancy (years) § Column 5: Mortality 2011 -10 -03 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 45

Exercises: Datasets – countries § Most often, we want to access complete columns in the matrix. This can be done by A(: , index) For example if you are interested in the lifeexpectancy column, it is recommended to do: >> life = x(: , 4); and then the vector life can be used to access the vector containing all life expectancies. 2011 -10 -03 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 46

Exercises: Datasets – countries § The sort() function can be used to sort all items of a vector in inclining order. >> life = A(: , 4); >> plot(life) 2011 -10 -03 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 47

Exercises: Datasets – countries § The sort() function can be used to sort all items of a vector in inclining order. >> life = A(: , 4); >> life. S = sort(life); >> plot(life. S) 2011 -10 -03 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 48

Exercises: Datasets – countries § The histogram hist() is useful for getting the distribution of the values of a vector. >> life = A(: , 4); >> hist(life) 2011 -10 -03 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 49

Exercises: Datasets – countries § Alternatively, a second parameter specifies the number of bars: >> life = A(: , 4); >> hist(life, 30) 2011 -10 -03 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 50

Exercise 1 § Statistics: Generate a vector of N random numbers with randn(N, 1) Calculate the mean and standard deviation. Do the mean and standard deviation converge to certain values, for an increasing N? Optional: Display the histogram and compare the output of the following two commands § hist(randn(N, 1)) § hist(rand(N, 1)) 2011 -10 -03 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 51

Exercise 2 § Demographics: From the countries. m dataset, find out why there is such a large difference between the mean and the median population of all countries. Hint: Use hist(x, n) Also sort() can be useful. § Plus: play with subplot() 2011 -10 -03 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 52

Exercise 3 § Demographics: From the countries. m dataset, see which columns have strongest correlation. Can you explain why these columns have stronger correlations? Hint: Use corrcoef() to find the correlation between columns. Use imagesc()to get an immediate visualization of the correlations. 2011 -10 -03 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 53

Exercise 4 – optional § Zipf’s law: Zipf’s law says that the rank, x, of cities (1: largest, 2: 2 nd largest, 3: 3 rd largest, . . . ) and the size, y, of cities (population) has a power-law relation: y ~ xb § Test if Zipf’s law holds for the three cases in the cities. m file. Try to estimate b. Hint: Use log() and plot() (or loglog()) § Plus: use the fitting tool cftool() to get b 2011 -10 -03 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 54

References § Additional material on plotting on the home page of the course http: //www. soms. ethz. ch/matlab § "Log or Linear? Distinct Intuitions of the Number Scale in Western and Amazonian Indigene Cultures”, Science 320 (5880): 1217 § http: //en. wikipedia. org/wiki/Weber%E 2%80%93 F echner_law 2011 -10 -03 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 55