Data Analytics for Sense-Making

Data Analytics for Sense-Making Vik Gopal, Lim Tiong Wee 2025-10-27 Contents Preface 6 1 Introduction to Python 7 1.1 Introduction...................................... 7 1.2 Installing Python and Jupyter Lab . . . . . . . . . . . . . . . . . . . . . . . . . 7 1.3 Basic Data Structures in Python . . . . . . . . . . . . . . . . . . . . . . . . . . 9 1.4 Slice Operator in Python . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 1.5 LoopsinPython ................................... 10 1.6 Strings......................................... 11 1.7 Functions, Modules and Packages . . . . . . . . . . . . . . . . . . . . . . . . . . 13 1.8 Object-Oriented Programming . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 1.9 Numpy......................................... 14 1.9.1 ArrayCreation................................ 15 1.9.2 Slice Operator in Multiple Dimensions . . . . . . . . . . . . . . . . . . . 16 1.9.3 BasicOperations............................... 17 1.9.4 Axis-wise Operations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 1.10Pandas......................................... 19 1.10.1Series ..................................... 19 1.10.2DataFrames.................................. 21 1.10.3ReadinginData ............................... 21 1.10.4BasicSelection ................................ 22 1.10.5 Indexing and Selecting Data . . . . . . . . . . . . . . . . . . . . . . . . . 22 1.10.6FilteringData ................................ 24 1.10.7MissingValues ................................ 25 1.11References....................................... 25 2 Statistical Inference 27 2.1 Introduction...................................... 27 2.1.1 HypothesisTests............................... 27 2.1.2 Confidence Intervals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 2.2 ComparingMeans .................................. 29 2.2.1 2-sampleTests ................................ 29 2.3 PairedSampleTests ................................. 33 2.3.1 FormalSet-up ................................ 34 2.4 ANoVA ........................................ 36 2.4.1 FormalSet-up ................................ 38 2.4.2 𝐹-TestinOne-WayANOVA ........................ 39 2.4.3 Assumptions ................................. 39 2.4.4 Comparing specific groups . . . . . . . . . . . . . . . . . . . . . . . . . . 41 2.4.5 Contrast Estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 2.4.6 Multiple Comparisons . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 2.5 CategoricalVariables................................. 44 2.5.1 𝜒2-Test for Independence . . . . . . . . . . . . . . . . . . . . . . . . . . 45 2 2.5.2 Measures of Association . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 2.5.3 OddsRatio .................................. 48 2.5.4 For Ordinal Variables . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 2.6 Summary ....................................... 51 2.7 References....................................... 51 2.7.1 WebsiteReferences.............................. 51 2.7.2 Documentation links . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 3 Unsupervised Learning 53 3.1 Introduction...................................... 53 3.2 Principal Components Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 3.2.1 FormalSet-up ................................ 55 3.3 Clustering....................................... 59 3.3.1 Hierarchical Clustering . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 3.3.2 Determining the optimal number of clusters . . . . . . . . . . . . . . . . 63 3.4 OutlierDetection................................... 65 3.5 Visualisation ..................................... 67 3.5.1 MDS...................................... 67 3.5.2 t-SNE ..................................... 69 3.6 References....................................... 72 3.6.1 Websitereferences .............................. 72 3.6.2 Videoreferences ............................... 72 4 Natural Language Processing 73 4.1 Introduction...................................... 73 4.2 Definitions....................................... 74 4.3 Overview of Applications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 4.4 TextPre-processing ................................. 75 4.4.1 Pre-processing Text with Gensim . . . . . . . . . . . . . . . . . . . . . . 76 4.5 RepresentationofText................................ 78 4.5.1 Sparse embeddings with Tf-idf . . . . . . . . . . . . . . . . . . . . . . . 78 4.5.2 Cosinesimilarity ............................... 81 4.5.3 DenseEmbeddings.............................. 82 4.6 Visualisation with t-SNE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 4.7 NeuralLanguageModels............................... 85 4.8 Applications...................................... 87 4.8.1 SentimentAnalysis.............................. 88 4.8.2 Information Retrieval . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 4.8.3 TopicModeling................................ 92 4.9 Interpretation of Neural Models . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 4.10References....................................... 95 4.10.1Videoexplainers ............................... 95 4.10.2WebsiteReferences.............................. 97 5 Linear Regression 98 5.1 Introduction...................................... 98 5.2 Simple Linear Regression . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 5.2.1 FormalSet-up ................................ 101 5.2.2 Estimation .................................. 102 5.2.3 Hypothesis Test for Model Significance . . . . . . . . . . . . . . . . . . . 103 5.2.4 Coefficient of Determination, 𝑅2...................... 104 3 5.3 Multiple Linear Regression . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106 5.3.1 FormalSetup................................. 106 5.3.2 Estimation .................................. 107 5.3.3 Adjusted 𝑅2................................. 107 5.3.4 HypothesisTests............................... 107 5.4 Including a Categorical Variable . . . . . . . . . . . . . . . . . . . . . . . . . . 111 5.4.1 Including an Interaction Term . . . . . . . . . . . . . . . . . . . . . . . . 113 5.5 ResidualAnalysis................................... 114 5.5.1 Standardised Residuals . . . . . . . . . . . . . . . . . . . . . . . . . . . 114 5.5.2 Scatterplots.................................. 116 5.5.3 InfluentialPoints............................... 117 5.6 Transformation.................................... 119 5.7 Summary, Further topics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120 5.8 References....................................... 120 5.8.1 WebsiteReferences.............................. 120 6 Time Series Analysis 121 6.1 Exploring Time Series Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121 6.2 Decomposing Time Series Data . . . . . . . . . . . . . . . . . . . . . . . . . . . 126 6.3 Forecasting ...................................... 130 6.3.1 Benchmarkmethods............................. 130 6.3.2 ARIMAModels................................ 133 6.3.3 ETS...................................... 138 6.3.4 ThetaMethod ................................ 142 6.3.5 Forecasting with Seasonal Decomposition . . . . . . . . . . . . . . . . . 145 6.4 MiscellaneousTopics................................. 146 6.4.1 Time Series Clustering . . . . . . . . . . . . . . . . . . . . . . . . . . . . 146 6.5 Summary ....................................... 150 6.6 References....................................... 150 6.6.1 Statsmodelspages.............................. 150 6.6.2 Forecasting principles and practice . . . . . . . . . . . . . . . . . . . . . 151 7 Simulation 152 7.1 RandomVariables .................................. 152 7.1.1 Discrete Random Variables . . . . . . . . . . . . . . . . . . . . . . . . . 152 7.1.2 Continuous Random Variables . . . . . . . . . . . . . . . . . . . . . . . 154 7.1.3 Generating Random Variates . . . . . . . . . . . . . . . . . . . . . . . . 155 7.2 General Principles in Simulation Studies . . . . . . . . . . . . . . . . . . . . . . 156 7.2.1 Introduction ................................. 156 7.2.2 Steps in a Simulation Study . . . . . . . . . . . . . . . . . . . . . . . . . 157 7.2.3 Theory .................................... 157 7.3 Object-Oriented Programming in Python . . . . . . . . . . . . . . . . . . . . . 158 7.4 Introduction to Agent Based Models . . . . . . . . . . . . . . . . . . . . . . . . 159 7.4.1 Introduction to Mesa . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 160 7.5 Agent and Model Classes (v1) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 160 7.5.1 DataCollection................................ 162 7.6 MultipleIterations .................................. 164 7.7 Agent and Model Classes (v2) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 165 7.7.1 Assessment of Income Inequality . . . . . . . . . . . . . . . . . . . . . . 165 7.7.2 DataCollection................................ 167 4 7.8 Agent and Model Classes (v3) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 168 7.8.1 Adding a spatial component . . . . . . . . . . . . . . . . . . . . . . . . . 168 7.8.2 DataCollection................................ 170 7.8.3 Boltzmann Model Final Comments . . . . . . . . . . . . . . . . . . . . . 172 7.9 Simpler Simulation Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 173 7.9.1 N-grammodels................................ 173 7.9.2 PowerAnalysis................................ 174 7.10 Summary: Simulation-Based Modeling . . . . . . . . . . . . . . . . . . . . . . . 176 7.11References....................................... 176 7.11.1MesaLinks .................................. 176 7.11.2 Other ABM Software . . . . . . . . . . . . . . . . . . . . . . . . . . . . 176 7.11.3 Reference Papers and Websites . . . . . . . . . . . . . . . . . . . . . . . 176 7.11.4 Other Simulation Software . . . . . . . . . . . . . . . . . . . . . . . . . 177 7.11.5OtherLinks.................................. 177 8 Supervised Learning 178 8.1 Introduction...................................... 178 8.2 Classification versus Regression . . . . . . . . . . . . . . . . . . . . . . . . . . . 178 8.3 Supervised Learning Workflow . . . . . . . . . . . . . . . . . . . . . . . . . . . 178 8.4 Scikit-learn ...................................... 179 8.4.1 Input Data Structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . 180 8.5 Measures of Performance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 181 8.5.1 ForClassification............................... 181 8.5.2 ForRegression ................................ 182 8.6 Classification ..................................... 183 8.6.1 DecisionTree................................. 184 8.6.2 Variable importance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 188 8.6.3 RandomForest................................ 191 8.7 Regression....................................... 195 8.7.1 Random Forest Regressor . . . . . . . . . . . . . . . . . . . . . . . . . . 196 8.8 Interpretability of Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 197 8.8.1 LIME ..................................... 197 8.8.2 ICEplots ................................... 200 8.9 Summary ....................................... 200 8.10References....................................... 201 8.10.1 Website and video references . . . . . . . . . . . . . . . . . . . . . . . . 201 8.10.2 Documentation references . . . . . . . . . . . . . . . . . . . . . . . . . . 201 9 Computer Vision 202 9.1 Introduction...................................... 202 9.2 ImageProcessing................................... 202 9.2.1 ReadingImages................................ 202 9.3 WorkingwithMasks................................. 203 9.4 Modifying Perspective of Images . . . . . . . . . . . . . . . . . . . . . . . . . . 206 9.5 ComputerVisionTasks ............................... 209 9.6 References....................................... 210 9.6.1 Opencv documentation . . . . . . . . . . . . . . . . . . . . . . . . . . . 210 9.6.2 Books ..................................... 210 9.6.3 Github repositories . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 210 Academic References 211 5 Preface IND5003 is a foundation course for students in a Master’s program for working adults, at the National University of Singapore. The title of the course is Data Analytics for Sense-Making The content of the course aims to bring students up to speed in a variety of analytic techniques that will be useful for them in later courses of the program. You can read some reviews of the course on NUSmods. In designing this course, we tried to cover the techniques from a couple of different angles: • From the viewpoint of the types of data encountered, e.g. time series, unstructured text data, image data, and so on. • From the viewpoint of the types of questions that could be asked of the data, e.g. are we attempting to plan for unseen scenarios (simulation), are we interrogating the data for hidden structure (unsupervised learning), and so on. The entire course uses the Python programming language. By the end of the course, the aim was for students to become familiar with Python for data analysis. If you are a fellow instructor and you find something useful in this textbook, please do let me know at [email protected]. If you need more details about anything, do feel free to write as well. So long, and thanks for reading! Vik https://blog.nus.edu.sg/stavg © 2025 Vik Gopal, Tiong Wee Lim. This book is licensed under Creative Commons Attribution 4.0 International (CC BY-SA 4.0). You are free to share and adapt with attribution. Code snippets are licensed under the MIT License. See the repository’s LICENSE and LICENSECODE. 6 1 Introduction to Python 1.1 Introduction Python is a general-purpose programming language. It is a higher-level language than C, C++ and Java in the sense that a Python program does not have to be compiled before execution. It was originally conceived back in the 1980s by Guido van Rossum at Centrum Wiskunde & Informatica (CWI) in the Netherlands. The language is named after a BBC TV show (Guido’s favorite program) “Monty Python’s Flying Circus”. Python reached version 1.0 in January 1994. Python 2.0 was released on October 16, 2000. Python 3.0, which is backwards-incompatible with earlier versions, was released on 3 December 2008. Python is a very flexible language; it is simple to learn yet is fast enough to be used in production. Over the past ten years, more and more comprehensive data science toolkits (e.g. scikit-learn, NTLK, tensorflow, keras) have been written in Python and are now the standard frameworks for those models. Python is an open-source software and thus it is free to use and extend. 1.2 Installing Python and Jupyter Lab To install Python, navigate to the official Python download page to obtain the appropriate installer for your operating system. ĺImportant For our class, please ensure that you are using Python 3.10.12. The next step is to create a virtual environment for this course. Virtual environments are specific to Python. They allow you to retain multiple versions of Python, and of packages, on the same computer. Go through the videos on Canvas relevant to your operating system to create a virtual environment and install Jupyter Lab on your machine. Jupyter notebooks are great for interactive work with Python, but more advanced users may prefer a full-fledged IDE. If you are an advanced user, and are comfortable with an IDE of your own choice (e.g. Spyder or VSCode), feel free to continue using that to run the codes for this course. ĺImportant Even if you are using Anaconda/Spyder/VSCode, you still need to create a virtual environment. 7 Jupyter notebooks consist of cells, which can be of three main types: • code cells, • output cells, and • markdown cells. Figure 1.1: Jupyter Lab In Figure 1.1, the red box labelled 1 is a markdown cell. It can be used to contain descriptions or summary of the code. The cells in the box labelled 2 are code cells. To run the codes from our notes, you can copy and paste the codes into a new cell, and then execute them with Ctrl-Enter. Try out this Easter egg that comes with any Python installation: import this The Zen of Python, by Tim Peters Beautiful is better than ugly. Explicit is better than implicit. Simple is better than complex. Complex is better than complicated. Flat is better than nested. Sparse is better than dense. Readability counts. Special cases aren't special enough to break the rules. 8 Although practicality beats purity. Errors should never pass silently. Unless explicitly silenced. In the face of ambiguity, refuse the temptation to guess. There should be one-- and preferably only one --obvious way to do it. Although that way may not be obvious at first unless you're Dutch. Now is better than never. Although never is often better than *right* now. If the implementation is hard to explain, it's a bad idea. If the implementation is easy to explain, it may be a good idea. Namespaces are one honking great idea -- let's do more of those! More information on using Jupyter notebooks can be obtained from this website. 1.3 Basic Data Structures in Python The main data structures in native1Python are •Lists, which are defined with [ ]. Lists are mutable. •Tuples, which are defined with ( ). Tuples are immutable. •Dictionaries, which are defined with { }. Dictionaries have keys and items. They are also mutable. Very soon, we shall see that for data analysis, the more common objects we shall deal with are dataframes (from pandas) and arrays (from numpy). However, the latter two require add-on packages; the three object classes listed above are baked into Python. Here is how we create lists, tuples and dictionaries. x_list =[1,2,3] x_tuple =(1,2,3) x_dict ={'a':1,'b':2,'c':3}# access with x_dict['a'] By the way, this is what mean by (im)mutable. The following is OK, because “x” is a list, and hence mutable: x=[1,3,5,7,8,9,10] x[3]=17 print(x) [1, 3, 5, 17, 8, 9, 10] But the following code will return an error, because x_tuple is a tuple, and hence immutable. x_tuple =(1,3,5,6,8,9,10) x_tuple[3]=17 1i.e., Python without any packages imported. 9 arr_real =np.linspace(start =0.2, stop =3.3, num =24).reshape((2,3,4)) arr_real array([[[0.2 , 0.33478261, 0.46956522, 0.60434783], [0.73913043, 0.87391304, 1.00869565, 1.14347826], [1.27826087, 1.41304348, 1.54782609, 1.6826087 ]], [[1.8173913 , 1.95217391, 2.08695652, 2.22173913], [2.35652174, 2.49130435, 2.62608696, 2.76086957], [2.89565217, 3.03043478, 3.16521739, 3.3 ]]]) Sometimes we need to create a placeholder array with the appropriate dimensions, and then fill it in later. This is preferable to growing an array by appending to it. np.zeros((3,5)) array([[0., 0., 0., 0., 0.], [0., 0., 0., 0., 0.], [0., 0., 0., 0., 0.]]) Similarly, there is also an np.ones() function. Instead of specifying the dimensions of an array ourselves, we can create arrays of zeros or ones in the shape of other existing arrays. The following code creates an array of ones, of the same shape as arr_real. np.ones_like(arr_real) array([[[1., 1., 1., 1.], [1., 1., 1., 1.], [1., 1., 1., 1.]], [[1., 1., 1., 1.], [1., 1., 1., 1.], [1., 1., 1., 1.]]]) 1.9.2 Slice Operator in Multiple Dimensions Multidimensional NumPy arrays can be accessed with comma separated slice notation. When fewer indices are provided than the number of axes, the missing indices are considered complete slices for the remaining dimensions. By the way, when printing, the last axis will be printed left-to-right, and the second last axis will be printed from top-to-bottom. The remaining axes will be printed with a line in between: arr_real 16 array([[[0.2 , 0.33478261, 0.46956522, 0.60434783], [0.73913043, 0.87391304, 1.00869565, 1.14347826], [1.27826087, 1.41304348, 1.54782609, 1.6826087 ]], [[1.8173913 , 1.95217391, 2.08695652, 2.22173913], [2.35652174, 2.49130435, 2.62608696, 2.76086957], [2.89565217, 3.03043478, 3.16521739, 3.3 ]]]) Here are examples based on this array. Try to guess what each will return before you run it: arr_real[1,2,3] arr_real[0,2, ::-1] arr_real[1,0:3:2] arr_real[:,2,:] An alternative to using integers in the selection index is to use Boolean indexing. This means that we use an array of True and False entries to determine which elements to return. arr_real >3 array([[[False, False, False, False], [False, False, False, False], [False, False, False, False]], [[False, False, False, False], [False, False, False, False], [False, True, True, True]]]) arr_real[arr_real >3] array([3.03043478, 3.16521739, 3.3 ]) 1.9.3 Basic Operations The next few commands generate 𝑈𝑛𝑖𝑓(0,1)random variates and store them in rectangular arrays. We shall see more about these particular functions in the chapter on Simulation (Section 7.1), but for now, take note that the number 1361 is known as a seed; it ensures the reproducibility of random numbers across sessions. # Setting a seed allows for reproducibility of random number generation # across sessions. rng =np.random.default_rng(1361) a=rng.uniform(size=(3,5)) b=rng.uniform(size=(3,5)) b array([[0.47031767, 0.43520574, 0.74929584, 0.19794778, 0.91397998], [0.95979996, 0.10301973, 0.82972687, 0.61118187, 0.04727632], [0.51642555, 0.32188431, 0.57800994, 0.73144202, 0.74020866]]) 17 The +operator, applied to conformable arrays, performs element-wise addition. a+b array([[0.92582851, 0.90610288, 1.3360719 , 0.43528561, 1.30890054], [1.78346478, 0.31376731, 1.51889626, 1.40816571, 0.21158293], [0.52725269, 0.57843605, 1.07039529, 1.54221413, 0.82312206]]) Similarly, the *operator performs element-wise multiplication. a*b array([[0.2142348 , 0.20493714, 0.43966886, 0.0469805 , 0.36094949], [0.79055346, 0.02171116, 0.57182236, 0.48710208, 0.00776781], [0.00559141, 0.08257998, 0.28460362, 0.59303279, 0.06137322]]) For matrix multiplication, we have to transpose b. a.dot(b.T) array([[1.26677078, 1.13630182, 1.19189672], [1.30342857, 1.87895687, 1.5961133 ], [0.721959 , 0.94481619, 1.02718103]]) 1.9.4 Axis-wise Operations Numpy and Pandas have several facilities for applying functions across array axes. It is important to be aware of how they work. For instance, this takes the mean across the 0th (first) axis. The resulting array has shape (3,4). arr_real.mean(axis =0) array([[1.00869565, 1.14347826, 1.27826087, 1.41304348], [1.54782609, 1.6826087 , 1.8173913 , 1.95217391], [2.08695652, 2.22173913, 2.35652174, 2.49130435]]) The top-left element comes from the average of arr_real[0,0,0] and arr_real[1,0,0]. Similarly, the element to the right of it comes from the average of arr_real[0,0,1] and arr_real[1,0,1]: (arr_real[0,0,1]+arr_real[1,0,1]) /2 1.143478260869565 Note that arr_real[0] is a 2D array, with shape (3, 4). Suppose we wish to compute the row means. This means we have to apply the operation by the column axis (axis = 1). 18 arr_real[0].mean(axis =1) array([0.40217391, 0.94130435, 1.48043478]) Here is a table with some common operations that you can apply on a numpy array. Method Description shape Returns dimensions, e.g. matrix1.shape TTransposes the array, e.g. matrix1.T mean Computes color row-wise means, e.g. matrix1.mean(axis=0) or matrix1.mean(axis=1) sum Computes color row-wise means, e.g. matrix1.sum(axis=0) or matrix1.sum(axis=1) argmax Return the index corresponding to the max within the specified dimension, e.g. matrix1.argmax(axis=0) for the position with the max within each column. reshape To change the dimensions, e.g. array1.reshape((5,1)) converts the array into a 5x1 matrix For instance, if we wanted to identify the row with the largest mean, we use argmax() on the resulting array. arr_real[0].mean(axis=1).argmax() 2 1.10 Pandas 1.10.1 Series ASeries is a one-dimensional labeled array. The axis labels are referred to as the index. The simplest way to create a Series is to pass a sequence and an index to pd.Series(). Example 1.1 (Creating Pandas Series).Consider the following data, from the football league in Spain. year =pd.Series(list(range(2010,2013) ) *3) team =["Barcelona","RealMadrid","Valencia"]*3 team.sort() team =pd.Series(team) wins =pd.Series([30,28,32,29,32,26,21,17,19]) draws =pd.Series([6,7,4,5,4,7,8,10,8]) losses =pd.Series([2,3,2,4,2,5,9,11,11]) #wins.index #wins.values 19 To access particular values, we can use the slice operator. wins[0:6:2] 0 30 2 32 4 32 dtype: int64 To convert a Series object to an ndarray, we use the following method: wins.to_numpy() array([30, 28, 32, 29, 32, 26, 21, 17, 19], dtype=int64) If we specify an index, we can use it to access values in the Series. s=pd.Series(rng.uniform(size=5), index=['a','b','c','d','e']) # s # s.index # s.values s[['a','c']] a 0.974173 c 0.948907 dtype: float64 Be careful when you combine the slice operator with label-based indexing. Unlike vanilla Python, Pandas includes both end-points! s['a':'d'] a 0.974173 b 0.270230 c 0.948907 d 0.675365 dtype: float64 20 1.10.2 DataFrames ADataFrame is a 2-dimensional labeled data structure with possibly different data types. It is the most commonly used Pandas object. The index of a DataFrame refers to the row labels (axis 0). The columns refer to the column labels (axis 1). DataFrames can be constructed from Series, dictionaries, lists and 2-d arrays. For our course, we will typically create a DataFrame directly from a file. Example 1.2 (Creating Pandas DataFrame from Series).We can create a DataFrame from the earlier series. laliga =pd.DataFrame({'Year': year, 'Team': team, 'Wins': wins, 'Draws': draws, 'Losses': losses }) To inspect a DataFrame, we can use info(),head() and tail() methods. laliga.head() Year Team Wins Draws Losses 0 2010 Barcelona 30 6 2 1 2011 Barcelona 28 7 3 2 2012 Barcelona 32 4 2 3 2010 RealMadrid 29 5 4 4 2011 RealMadrid 32 4 2 1.10.3 Reading in Data Pandas can read in data stored in multiple formats, including CSV, tab-separated files, and Excel files. Example 1.3 (Happiness dataset).The CSV file read in here contains the happiness scores of 164 countries from 2015 to 2017. Click here for a full report on the dataset. The final score was based on many other factors (such as GDP per capita, family, freedom etc) which is included in the file as well. In each year, not all of the 164 countries had their scores surveyed and taken. This results in some countries having missing values (NaN) in certain years. happ =pd.read_csv('data/happiness_report.csv', header=0, na_values='NA') 21 1.10.4 Basic Selection In dataframes, row selection can technically be done with integers with the slice operator. In practice, this is not used, because we typically wish to select a set of rows based on a condition. print(happ[10:12]) Country Happiness.Rank Happiness.Score GDP Family \ 10 Israel 11.0 7.278 1.22857 1.22393 11 Costa Rica 12.0 7.226 0.95578 1.23788 Life.Expectancy Freedom Govt.Corruption Generosity Dystopia.Residual \ 10 0.91387 0.41319 0.07785 0.33172 3.08854 11 0.86027 0.63376 0.10583 0.25497 3.17728 Year 10 2015 11 2015 To select columns, you may use a list of column names. happ[['GDP','Freedom']].head() # note the difference with happ['GDP'] GDP Freedom 0 1.39651 0.66557 1 1.30232 0.62877 2 1.32548 0.64938 3 1.45900 0.66973 4 1.32629 0.63297 ÁWarning Remember that we are not working with numpy arrays, so this will not work: happ[0:10,2:4] 1.10.5 Indexing and Selecting Data The two main methods of advanced data selection use the .loc and .iloc functions. Although we call them functions, they are summoned using the [ ] notation. The .loc is primarily label-based. The common allowed inputs to .loc are • a single label, • a list of labels, • a slice object, 22 • a boolean array. The .iloc is primarily an integer-based input. The common allowed inputs to .iloc are • a single integer, • a list of integers, • a slice object, • a boolean array. When selecting from a DataFrame with .loc or .iloc, we can provide a comma-separated index, just as with NumPy. It is good to keep this reference page bookmarked. Take note that this next command will only work if the index is made up of integers! print(happ.loc[2:5]) Country Happiness.Rank Happiness.Score GDP Family \ 2 Denmark 3.0 7.527 1.32548 1.36058 3 Norway 4.0 7.522 1.45900 1.33095 4 Canada 5.0 7.427 1.32629 1.32261 5 Finland 6.0 7.406 1.29025 1.31826 Life.Expectancy Freedom Govt.Corruption Generosity Dystopia.Residual \ 2 0.87464 0.64938 0.48357 0.34139 2.49204 3 0.88521 0.66973 0.36503 0.34699 2.46531 4 0.90563 0.63297 0.32957 0.45811 2.45176 5 0.88911 0.64169 0.41372 0.23351 2.61955 Year 2 2015 3 2015 4 2015 5 2015 Notice below how the slice operator is inclusive when we use .loc, but not inclusive when we use .iloc. Consider the following output, using .loc. happ.loc[2:10:4,"GDP":"Generosity":2] GDP Life.Expectancy Govt.Corruption 2 1.32548 0.87464 0.48357 6 1.32944 0.89284 0.31814 10 1.22857 0.91387 0.07785 In order to obtain the same output with .iloc, we have to use this instead: happ.iloc[2:11:4,3:8:2] 23 GDP Life.Expectancy Govt.Corruption 2 1.32548 0.87464 0.48357 6 1.32944 0.89284 0.31814 10 1.22857 0.91387 0.07785 1.10.6 Filtering Data Suppose we are interested in the very happy countries. Here is how we can filter the data with a boolean array. Example 1.4 (Happiest countries).Suppose we retrieve the happiest countries; those with a score more than 6.95. Can you surmise why we use this value? happiest =happ[happ['Happiness.Score']>6.95] happiest.Country.unique() array(['Switzerland', 'Iceland', 'Denmark', 'Norway', 'Canada', 'Finland', 'Netherlands', 'Sweden', 'New Zealand', 'Australia', 'Israel', 'Costa Rica', 'Austria', 'Mexico', 'United States', 'Brazil', 'Ireland', 'Germany'], dtype=object) Notice that there isn’t a single Asian or African country in the happiest 10% of countries! When filtering, we can also combine Boolean indices. # Top 3 happiest countries in 2015 print(happ[(happ.Year == 2015)&(happ['Happiness.Rank']<= 3)]) Country Happiness.Rank Happiness.Score GDP Family \ 0 Switzerland 1.0 7.587 1.39651 1.34951 1 Iceland 2.0 7.561 1.30232 1.40223 2 Denmark 3.0 7.527 1.32548 1.36058 Life.Expectancy Freedom Govt.Corruption Generosity Dystopia.Residual \ 0 0.94143 0.66557 0.41978 0.29678 2.51738 1 0.94784 0.62877 0.14145 0.43630 2.70201 2 0.87464 0.64938 0.48357 0.34139 2.49204 Year 0 2015 1 2015 2 2015 24 1.10.7 Missing Values The .info() method will yield information on missing values, column by column. We can see there are 21 rows with missing values. happ.info() <class 'pandas.core.frame.DataFrame'> RangeIndex: 492 entries, 0 to 491 Data columns (total 11 columns): # Column Non-Null Count Dtype --- ------ -------------- ----- 0 Country 492 non-null object 1 Happiness.Rank 471 non-null float64 2 Happiness.Score 471 non-null float64 3 GDP 471 non-null float64 4 Family 471 non-null float64 5 Life.Expectancy 471 non-null float64 6 Freedom 471 non-null float64 7 Govt.Corruption 471 non-null float64 8 Generosity 471 non-null float64 9 Dystopia.Residual 471 non-null float64 10 Year 492 non-null int64 dtypes: float64(9), int64(1), object(1) memory usage: 42.4+ KB Sometimes, it is appropriate to drop rows with missing values. This can be done with the .dropna method. Remember that it returns a new dataframe. The original one remains unchanged, unless you include the inplace=True argument. new_df =happ.dropna() # pd.isna(happ) 1.11 References In this chapter, we have introduced the following data science tools: • Python programming language • Jupyter notebooks for working with Python • Computer set-up for data science with Python Python is very widely used for data science, and especially for the machine learning aspect of it. For those of you with intentions to take up the GC in Deep Learning or Data Mining, it will be critical to be familiar with the language. It will be used again in at least DSA5102. 1. Regular expression HOWTO A tutorial with examples on regular expressions (for manipulating and searching through strings). 2. Formatting string literals: Or just f-strings 25 A Quantile-Quantile plot is a graphical diagnostic tool for assessing if a dataset follows a particular distribution. Most of the time we would be interested in comparing against a Normal distribution. A QQ-plot plots the standardized sample quantiles against the theoretical quantiles of a N(0; 1) distribution. If they fall on a straight line, then we would say that there is evidence that the data came from a normal distribution. Especially for unimodal datasets, the points in the middle will fall close to the line. The value of a QQ-plot is in judging if the tails of the data are fatter or thinner than the tails of the Normal. (a) Thinner (a) Fatter Example 2.3 (Abalone measurements QQ-plots).If we compare the qq-plots from the data (below) with the reference (above), we can infer that for females, the left tail is thinner than a Normal - it abruptly cuts off. For males, both the left and the right tail are fatter than a Normal’s. f, axs =plt.subplots(1,2, figsize=(8,4)) tmp =plt.subplot(121) sm.qqplot(x, line="q", ax=tmp) tmp.set_title('Females') tmp =plt.subplot(122) sm.qqplot(y, line="q", ax=tmp) tmp.set_title('Males'); 32 21012 Theoretical Quantiles 0.1 0.2 0.3 0.4 0.5 0.6 Sample Quantiles Females 21012 Theoretical Quantiles 0.1 0.2 0.3 0.4 0.5 0.6 Sample Quantiles Males There are also numerous hypothesis tests for Normality. If you are keen on learning about them, take a look at the references Section 2.7. We also need to assess if the variances are equal. While there are many hypothesis tests specifically for assessing if variances are equal (e.g. Levene, Bartlett), in our class, I advocate a simple rule of thumb. If the larger s.d is more than twice the smaller one, than we should not use the equal variance form of the test. This rule of thumb is widely used in practice (see the references Section 2.7). abl.groupby('gender').describe() viscera count mean std min 25% 50% 75% max gender F 50.0 0.28241 0.108707 0.095 0.201250 0.275 0.365125 0.575 M 50.0 0.30220 0.108746 0.040 0.253125 0.310 0.348750 0.638 We would conclude that there is no significant difference between the mean viscera weight of males and females. 2.3 Paired Sample Tests The data in a paired sample test also arises from two groups, but the two groups are not independent. A very common scenario that gives rise to this test is when the same subject receives both treatments. His/her measurement under each treatment gives rise to a measurement in each group. However, the measurements are no longer independent. 33 Example 2.4 (Reaction time of drivers).Consider a study on 32 drivers sampled from a driving school. Each driver is put in a simulation of a driving situation, where a target flashes red and green at random periods. Whenever the driver sees red, he/she has to press a brake button. For each driver, the study is carried out twice - at one of the repetitions, the individual carries on a phone conversation while at the other, the driver listens to the radio. Each measurement falls under one of two groups - “phone” or “radio”, but the measurements for driver 𝑖are clearly related. Some people might just have a slower/faster baseline reaction time! This is a situation where a paired sample test is appropriate, not an independent sample test. 2.3.1 Formal Set-up Suppose that we observe 𝑋1,…,𝑋𝑛independent observations from group 1 and 𝑌1,…,𝑌𝑛 independent observations from group 2. However the pair (𝑋𝑖,𝑌𝑖)are correlated. Similar to the previous section, it is assumed that 𝑋𝑖∼ 𝑁(𝜇1,𝜎2 1),𝑖=1,…,𝑛 (2.7) 𝑌𝑗∼ 𝑁(𝜇2,𝜎2 2),𝑗=1,…,𝑛 (2.8) We let 𝐷𝑖=𝑋𝑖−𝑌𝑖for 𝑖=1,…,𝑛. It follows that 𝐷𝑖∼𝑁(𝜇1−𝜇2,𝜎2 1+𝜎2 2−2𝑐𝑜𝑣(𝑋𝑖,𝑌𝑖)) The null and alternative hypotheses are stated in terms of the distribution of 𝐷𝑖: 𝐻0∶ 𝜇𝐷=0 𝐻1∶ 𝜇𝐷≠0 The test statistic for this test is: 𝑇2= 𝐷−0 𝑠/√𝑛 where 𝑠2=∑𝑛 𝑖=1(𝐷𝑖− 𝐷)2 (𝑛−1) Under 𝐻0, the test statistic 𝑇2∼𝑡𝑛−1. When we use a software to apply the test above, it will typically also return a confidence interval, computed as  𝐷±𝑡𝑛−1,1−𝛼/2×𝑠/√𝑛 34 Example 2.5 (Example: Heart Rate Before/After Treadmill).The following dataset comes from a textbook. In a self-recorded experiment, an individual recorded his heart rate before using a treadmill (baseline) and 5 minutes after use, for 12 days in 2006. hr_df =pd.read_csv("data/health_promo_hr.csv") p_test_out =stats.ttest_rel(hr_df.baseline, hr_df.after5) print(f""" The $p$-value for the test is {p_test_out.pvalue:.2g}. The difference in means is {hr_df.baseline.mean() -hr_df.after5.mean():.3f}. """) The $p$-value for the test is 2.2e-10. The difference in means is -15.229. While we do not include them here, it is imperative to also make the checks for Normality. If you were to make them, you would realise that the sample size is rather small. It is difficult to make the case for Normality here. Here is a plot that is particularly useful in paired sample studies: ax1 =hr_df.plot(x='baseline', y='after5',kind='scatter', marker='o', edgecolor='blue', color='none') group_means =hr_df.loc[:, ['baseline','after5']].mean(axis=0) ax1.set_xlim(75,105) ax1.set_ylim(75,105) ax1.plot([75,105], [75,105], color="lightblue", linestyle="dashed"); ax1.scatter(group_means.iloc[0], group_means.iloc[1], marker='o', edgecolor='blue', color='none', s=100); ax1.set_title('Agreement of after5 and baseline'); 35 75 80 85 90 95 100 105 baseline 75 80 85 90 95 100 105 after5 Agreement of after5 and baseline Figure 2.4: Agreement plot In an agreement plot, we check to see if the scatter of points is close to the diagonal. If it is, we have evidence of little difference between the means. The largest circle in the plot corresponds to the mean readings, before and after the treadmill run. 2.4 ANoVA In this section, we introduce the one-way analysis of variance (ANOVA), which generalises the 𝑡-test methodology to more than 2 groups. Hypothesis tests in the ANOVA framework require the assumption of Normality. While the 𝐹-test in ANOVA provides a determination of whether or not the group means are different, in practice, we would always want to follow up with specific comparisons between groups as well. Example 2.6 (Heifers dataset).The following example was taken from Introduction to Statistical Data Analysis for Life Sciences. An experiment with dung from heifers was carried out in order to explore the influence of antibiotics on the decomposition of dung organic material. As part of the experiment, 36 heifers were randomly assigned into six groups. Note that a heifer is a young, female cow that has not had her first calf yet. Antibiotics of different types were added to the feed for heifers in five of the groups. The remaining group served as a control group. For each heifer, a bag of dung was dug into the soil, and after 8 weeks the amount of organic material was measured for each bag. Here is a boxplot of the data from each group, along with summary statistics in a table below. 36 heifers =pd.read_csv('data/antibio.csv') sns.boxplot(heifers, x='type', y='org',); Ivermect Alfacyp Enroflox Spiramyc Fenbenda Control type 2.4 2.5 2.6 2.7 2.8 2.9 3.0 3.1 org Compared to the control group, it does appear that the median organic weight of the dung from the other heifer groups is higher. The following table displays the mean, standard deviation, and count from each group: heifers.groupby('type').describe() org count mean std min 25% 50% 75% max type Alfacyp 6.0 2.895000 0.116748 2.75 2.7950 2.915 2.9900 3.02 Control 6.0 2.603333 0.118771 2.43 2.5450 2.595 2.6825 2.76 Enroflox 6.0 2.710000 0.161988 2.42 2.6775 2.735 2.8075 2.88 Fenbenda 6.0 2.833333 0.123558 2.66 2.7675 2.850 2.8725 3.02 Ivermect 6.0 3.001667 0.109438 2.81 2.9625 3.045 3.0600 3.11 Spiramyc 4.0 2.855000 0.054467 2.80 2.8300 2.845 2.8700 2.93 Observe that the Spiramycin group only yielded 4 readings instead of 6. Our goal in this topic is to understand a technique for assessing if group means are statistically different from one another. Here are the specific analyses that we shall carry out: Heifers: Questions of Interest 1. Is there any significant difference, at 5% level, between the mean decomposition level of the groups? 2. At 5% level, is the mean level for Enrofloxacin different from the control group? 37 3. Pharmacologically speaking, Ivermectin and Fenbendazole are similar to each other. Let us call this sub-group (A). They work differently than Enrofloxacin. At 5% level, is there a significant difference between the mean from sub-group A and Enrofloxacin? 2.4.1 Formal Set-up Suppose there are 𝑘groups with 𝑛𝑖observations in the 𝑖-th group. The 𝑗-th observation in the 𝑖-th group will be denoted by 𝑌𝑖𝑗. In the One-Way ANOVA, we assume the following model: 𝑌𝑖𝑗 =𝜇+𝛼𝑖+𝑒𝑖𝑗,𝑖=1,…,𝑘,𝑗=1,…,𝑛𝑖(2.9) •𝜇is a constant, representing the underlying mean of all groups taken together. •𝛼𝑖is a constant specific to the 𝑖-th group. It represents the difference between the mean of the 𝑖-th group and the overall mean. •𝑒𝑖𝑗 represents random error about the mean 𝜇+𝛼𝑖for an individual observation from the 𝑖-th group. In terms of distributions, we assume that the 𝑒𝑖𝑗are i.i.d from a Normal distribution with mean 0 and variance 𝜎2. This leads to the model for each observation: 𝑌𝑖𝑗 ∼𝑁(𝜇+𝛼𝑖,𝜎2) It is not possible to estimate both 𝜇and all the 𝑘different 𝛼𝑖’s, since we only have 𝑘observed mean values for the 𝑘groups. For identifiability purposes, we need to constrain the parameters. There are two common constraints used, and note that different software have different defaults: 1. Setting ∑𝑘 𝑖=1𝛼𝑖=0, or 2. Setting 𝛼1=0. Continuing on from the equation for 𝑌𝑖𝑗, let us denote the mean for the 𝑖-th group as 𝑌𝑖, and the overall mean of all observations as 𝑌. We can then write the deviation of an individual observation from the overall mean as: 𝑌𝑖𝑗−𝑌 =(𝑌𝑖𝑗−𝑌𝑖) ⏟⏟⏟⏟⏟ within +(𝑌𝑖−𝑌) ⏟ between The first term on the right of the above equation is an indication of within-group variability. The second term on the right is an indication of between-group variability. The intuition behind the ANOVA procedure is that if the between-group variability is large and the within-group variability is small, then we have evidence that the group means are different. If we square both sides of the above equation and sum over all observations, we arrive at the following equation; the essence of ANOVA: 𝑘 ∑ 𝑖=1 𝑛𝑖 ∑ 𝑗=1(𝑌𝑖𝑗−𝑌)2=𝑘 ∑ 𝑖=1 𝑛𝑖 ∑ 𝑗=1(𝑌𝑖𝑗−𝑌𝑖)2+𝑘 ∑ 𝑖=1 𝑛𝑖 ∑ 𝑗=1(𝑌𝑖−𝑌)2 38 The squared sums above are referred to as: 𝑆𝑆𝑇=𝑆𝑆𝑊+𝑆𝑆𝐵 •𝑆𝑆𝑇: Sum of Squares Total, •𝑆𝑆𝑊: Sum of Squares Within, and •𝑆𝑆𝐵: Sum of Squares Between. In addition the following definitions are important for understanding the ANOVA output: 1. The Between Mean Square: 𝑀𝑆𝐵=𝑆𝑆𝐵 𝑘−1 2. The Within Mean Square: 𝑀𝑆𝑊=𝑆𝑆𝑊 𝑛−𝑘 The mean squares are estimates of the variability between and within groups. The ratio of these quantities is the test statistic. 2.4.2 𝐹-Test in One-Way ANOVA The null and alternative hypotheses are: 𝐻0∶ 𝛼𝑖=0for all 𝑖 𝐻1∶ 𝛼𝑖≠0for at least one 𝑖 The test statistic is given by 𝐹= 𝑀𝑆𝐵 𝑀𝑆𝑊 Under 𝐻0, the statistic 𝐹follows an 𝐹distribution with 𝑘−1and 𝑛−𝑘degrees of freedom. 2.4.3 Assumptions These are the assumptions that will need to be validated. 1. The observations are independent of each other. This is usually a characteristic of the design of the experiment, and is not something we can always check from the data. 2. The errors are Normally distributed. Residuals can be calculated as follows: 𝑌𝑖𝑗−𝑌𝑖 The distribution of these residuals should be checked for Normality. 3. The variance within each group is the same. In ANOVA, the 𝑀𝑆𝑊is a pooled estimate (across the groups) that is used; in order for this to be valid, the variance within each group should be identical. As in the 2-sample situation, we shall avoid separate hypotheses tests and proceed with the rule-of-thumb that if the ratio of the largest to smallest standard deviation is less than 2, we can proceed with the analysis. Example 2.7 (Heifers ANOVA F-test).Let us now apply the F-test to the heifers dataset. 39 heifer_lm =ols('org ~ type', data=heifers).fit() anova_tab =sm.stats.anova_lm(heifer_lm, type=3,) anova_tab df sum_sq mean_sq F PR(>F) type 5.0 0.590824 0.118165 7.972558 0.00009 Residual 28.0 0.415000 0.014821 NaN NaN At the 5% significance level, we reject the null hypothesis to conclude that the group means are significantly different from one another. This answers the first question we set out to. To extract the estimated parameters, we can use the following code: heifer_lm.summary() Dep. Variable: org R-squared: 0.587 Model: OLS Adj. R-squared: 0.514 Method: Least Squares F-statistic: 7.973 Date: Fri, 10 Oct 2025 Prob (F-statistic): 8.95e-05 Time: 09:42:27 Log-Likelihood: 26.655 No. Observations: 34 AIC: -41.31 Df Residuals: 28 BIC: -32.15 Df Model: 5 Covariance Type: nonrobust coef std err t P>|t|[0.025 0.975] Intercept 2.8950 0.050 58.248 0.000 2.793 2.997 type[T.Control] -0.2917 0.070 -4.150 0.000 -0.436 -0.148 type[T.Enroflox] -0.1850 0.070 -2.632 0.014 -0.329 -0.041 type[T.Fenbenda] -0.0617 0.070 -0.877 0.388 -0.206 0.082 type[T.Ivermect] 0.1067 0.070 1.518 0.140 -0.037 0.251 type[T.Spiramyc] -0.0400 0.079 -0.509 0.615 -0.201 0.121 Omnibus: 2.172 Durbin-Watson: 2.146 Prob(Omnibus): 0.338 Jarque-Bera (JB): 1.704 Skew: -0.545 Prob(JB): 0.427 Kurtosis: 2.876 Cond. No. 6.71 Notes: [1] Standard Errors assume that the covariance matrix of the errors is correctly specified. When estimating, Python sets one of the 𝛼𝑖to be equal to 0. We can tell from the output that the constraint has been placed on the coefficient for Alfacyp (since it is missing). From the output, we can read off (the Intercept term) that the estimate for Alfacyp is precisely 2.895+0=2.895 To check the assumptions, we can use the following code: 40 f, axs =plt.subplots(1,2, figsize=(10,4)) tmp =plt.subplot(121) heifer_lm.resid.hist(); tmp =plt.subplot(122) sm.qqplot(heifer_lm.resid, line="q", ax=tmp); 0.3 0.2 0.1 0.0 0.1 0.2 0 1 2 3 4 5 6 1.5 1.0 0.5 0.0 0.5 1.0 1.5 Theoretical Quantiles 0.3 0.2 0.1 0.0 0.1 0.2 Sample Quantiles 2.4.4 Comparing specific groups The 𝐹-test in a One-Way ANOVA indicates if all means are equal, but does not provide further insight into which particular groups differ. If we had specified beforehand that we wished to test if two particular groups 𝑖1and 𝑖2had different means, we could do so with a t-test. Here are the details to compute a confidence interval in this case: 1. Compute the estimate of the difference between the two means: 𝑌𝑖1−𝑌𝑖2 2. Compute the standard error of the above estimator: √ √ √ ⎷𝑀𝑆𝑊(1 𝑛𝑖1+1 𝑛𝑖2) 3. Compute the 100(1−𝛼)confidence interval as: 𝑌𝑖1−𝑌𝑖2±𝑡𝑛−𝑘,𝛼/2×√ √ √ ⎷𝑀𝑆𝑊(1 𝑛𝑖1+1 𝑛𝑖2) ĹNote If you notice from the summary statistics output for each group, the rule-of-thumb regarding standard deviations has not been satisfied. The ratio of largest to smallest standard deviations is slightly more than 2. Hence we should not continue with ANOVA; the pooled 41 2.5.3 Odds Ratio The most generally applicable measure of association, for 2x2 tables with nominal variables, is the Odds Ratio (OR). Suppose we have 𝑋and 𝑌to be Bernoulli random variables with (population) success probabilities 𝑝1and 𝑝2. We define the odds of success for 𝑋to be 𝑝1 1−𝑝1 Similarly, the odds of success for random variable 𝑌is 𝑝2 1−𝑝2. In order to measure the strength of their association, we use the odds ratio: 𝑝1/(1−𝑝1) 𝑝2/(1−𝑝2) The odds ratio can take on any value from 0 to ∞. • A value of 1 indicates no association between 𝑋and 𝑌. If 𝑋and 𝑌were independent, this is what we would observe. • Deviations from 1 indicate stronger association between the variables. • Note that deviations from 1 are not symmetric. For a given pair of variables, an association of 0.25 or 4 is the same - it is just a matter of which variable we put in the numerator odds. Due to the above asymmetry, we often use the log-odds-ratio instead: log 𝑝1/(1−𝑝1) 𝑝2/(1−𝑝2) • Log-odds-ratios can take values from −∞to ∞. • A value of 0 indicates no association between 𝑋and 𝑌. • Deviations from 0 indicate stronger association between the variables, and deviations are now symmetric; a log-odds-ratio of -0.2 indicates the same strength as 0.2, just the opposite direction. To obtain a confidence interval for the odds-ratio, we work with the log-odds ratio and then exponentiate the resulting interval. Here are the steps: 1. The sample data in a 2x2 table can be labelled as 𝑛11,𝑛12,𝑛21,𝑛22. 2. The sample odds ratio is  𝑂𝑅=𝑛11×𝑛22 𝑛12×𝑛21 3. For a large sample size, it can be shown that log  𝑂𝑅follows a Normal distribution. Hence a 95% confidence interval can be obtained through log 𝑛11×𝑛22 𝑛12×𝑛21 ±𝑧0.025×𝐴𝑆𝐸(log  𝑂𝑅) where • the ASE (Asymptotic Standard Error) of the estimator is √1 𝑛11 +1 𝑛12 +1 𝑛21 +1 𝑛22 48 Example 2.15 (Chest pain and gender odds ratio).Let us compute the confidence interval for the odds ratio in the chest pain and gender example from earlier. chest_tab2 =sm.stats.Table2x2(chest_array) chest_tab2.summary() Estimate SE LCB UCB p-value Odds ratio 1.353 0.863 2.123 0.188 Log odds ratio 0.303 0.230 -0.148 0.753 0.188 Risk ratio 1.322 0.872 2.004 0.188 Log risk ratio 0.279 0.212 -0.137 0.695 0.188 2.5.4 For Ordinal Variables When both variables are ordinal, it is often useful to compute the strength (or lack) of any monotone trend association. It allows us to assess if As the level of 𝑋increases, responses on 𝑌tend to increase toward higher levels, or responses on 𝑌tend to decrease towards lower levels. For instance, perhaps job satisfaction tends to increase as income does. In this section, we shall discuss a measure for ordinal variables, analogous to Pearson’s correlation for quantitative variables, that describes the degree to which the relationship is monotone. It is based on the idea of a concordant or discordant pair of subjects. • A pair of subjects is concordant if the subject ranked higher on 𝑋also ranks higher on 𝑌. • A pair is discordant if the subject ranking higher on 𝑋ranks lower on 𝑌. • A pair is tied if the subjects have the same classification on 𝑋and/or 𝑌. If we let •𝐶: number of concordant pairs in a dataset, and •𝐷: number of discordant pairs in a dataset. Then if 𝐶is much larger than 𝐷, we would have reason to believe that there is a strong positive association between the two variables. Here are two measures of association based on 𝐶and 𝐷: 1. Goodman-Kruskal 𝛾is computed as 𝛾=𝐶−𝐷 𝐶+𝐷 2. Kendall 𝜏𝑏is 𝜏𝑏=𝐶−𝐷 𝐴 where 𝐴is a normalising constant that results in a measure that works better with ties, and is less sensitive than 𝛾to the cut-points defining the categories. 𝛾has the advantage that it is more easily interpretable. 49 For both measures, values close to 0 indicate a very weak trend, while values close to 1 (or -1) indicate a strong positive (negative) association. Example 2.16 (Job satisfaction by income).Consider the following table, obtained from Agresti (2012). The original data come from a nationwide survey conducted in the US in 1996. us_svy_tab =np.array([[1,3,10,6], [2,3,10,7], [1,6,14,12], [0,1,9,11]]) col_names =['V. Diss','L. Diss','M. Sat','V. Sat'] row_names =['<15K','15-25K','25-40K','>40K'] sns.heatmap(us_svy_tab, annot=True, square=True, fmt='', xticklabels=col_names, yticklabels=row_names, cmap='Reds', cbar=False, ); V. Diss L. Diss M. Sat V. Sat <15K15-25K25-40K>40K 1 3 10 6 2 3 10 7 16 14 12 019 11 In measuring the association between these two variables (job satisfaction and income), we are interested in answering questions such as: If individual A has higher income than individual B, is individual A more likely to be satisfied in his/her job? 𝛾and 𝜏𝑏are measures that quantify this association. For the function in scipy.stats that computes this association, we need the data in “long format”. Hence we unroll it manually before summoning the method. dim1 =us_svy_tab.shape x=[];y=[] for iin range(0, dim1[0]): for jin range(0, dim1[1]): for kin range(0, us_svy_tab[i,j]): 50 x.append(i) y.append(j) stats.kendalltau(x, y, variant='b') SignificanceResult(statistic=0.15235215134659688, pvalue=0.0860294855671433) The output shows that 𝜏𝑏=0.15is positive, and is borderline significant at 5% level. We can conclude that there is a weak association between job satisfaction and income. 2.6 Summary In this topic, the primary take-aways are the notions of hypothesis tests (HT) and confidence intervals (CI). Both of these approaches aim to uncover information about a population using a sample. I strongly advocate the choice of CI over HT, as the latter option only provides a binary decision. CIs, on the other hand, provide a range of values to provide a better understanding of the situation. Using a small set of applications, we have demonstrated these techniques. We covered a common scenario where a researcher may need to compare the means between several groups. We then moved onto another common situation, where one might be faced with a contingency table. A common query is: What should we do if the assumptions of the test are not fulfilled? In those circumstances, one possibility is to turn to nonparametric versions of the test. For instance, the Kruskal-Wallis test is a HT for comparing the means of several groups whose distributions are not Normal. In the section on categorical data, we introduced measures of association for categorical variables - odds ratio, and Kendall 𝜏. When we have contingency tables, a common choice has been to display barcharts, or heat maps and use those to make decisions on. However, I encourage you to use the measures of association instead. They are intuitive to understand, and can allow you to compare sub-groups or variation of association over time. With HT and with CI, there has been a tendency to manipulate the data or tests until the desired outcome has been reached. Do be watchful of this. Conducting multiple tests increases the false positive error rate, so please do avoid this as well. If you conduct multiple tests, adjust for multiple testing. Use statitsical inference techniques as a guide together with common sense and domain expertise. 2.7 References 2.7.1 Website References 1. Inference recap from Penn State: •Hypothesis testing recap •Confidence intervals recap 2. Tests for Normality More information on the Kolmogorov-Smirnov and Shapiro-Wilks Tests for Normality. 51 3. Overview of 𝑡-tests This page includes the rule of thumb about deciding when to use the equal variance test, and when to use the unequal variances version. 2.7.2 Documentation links 1. statsmodels ANOVA A more complete example on the application of ANOVA. Return to this when we complete the topic on regression. 2. Useful functions from scipy.stats: Under the “Hypothesis Tests” section of this page, you can find: • Kruskal Wallis - the nonparametric equivalent of ANOVA • Wilcoxon signed rank - the nonparametric equivalent of paired sample t-test • Mann Whitney test - the nonparametric equivalent of independent samples t test • Somer’s D - measure of association for two categorical variables (one ordinal and one nominal). 52 3 Unsupervised Learning 3.1 Introduction Suppose that we have a set of 𝑁observations (𝑥1,𝑥2,…,𝑥𝑁)of a random 𝑝-vector 𝑋. The goal in unsupervised learning is to infer properties of the probability density of 𝑋. Note the primary difference with supervised learning (Section 8.1). In that context, we will have a set of labels 𝑦1,…,𝑦𝑁in addition to the 𝑥𝑖’s. Here, we do not have labelled data. In situations where 𝑝≤3, then graphical methods and numerical summaries such as correlations will suffice to help us understand the structure of the data. However, these methods breakdown as soon as 𝑝increases beyond 3. This topic introduces techniques that we can use, even when 𝑝is large, to : 1. Understand and interpret the main sources of variation in the data, 2. Identify “groups” or clusters within the data for further study. 3. Visualise high-dimensional data In the code below, one of the modules imported is clust. It contains a function that we shall use in the section on clustering (Section 3.3). import pandas as pd import numpy as np from scipy.cluster import hierarchy import seaborn as sns import matplotlib.pyplot as plt import plotly.express as px from itables import show from pprint import pprint import folium import geopandas from sklearn import decomposition, preprocessing from sklearn.metrics import pairwise_distances from sklearn.manifold import MDS, TSNE from sklearn.ensemble import IsolationForest from sentence_transformers import SentenceTransformer from ind5003 import clust 53 Example 3.1 (Wine quality data).The UCI Machine Learning Repository contains a dataset on Wine Quality. It consists of two tables - one corresponding to white wine and one corresponding to red wine. Each table contains the following columns: 1. fixed acidity 2. volatile acidity 3. citric acid 4. residual sugar 5. chlorides 6. free sulfur dioxide 7. total sulfur dioxide 8. density 9. pH 10. sulphates 11. alcohol 12. quality (score between 0 and 10) Columns 1 - 11 are numeric variables, measured objectively on the wines. Column 12 is a subjective evaluation made by wine experts, based on sensory data. Each quality score is the median of at least 3 evaluations. Although this dataset was created for a supervised learning problem, we shall use it to practice unsupervised learning techniques. To do so, we shall ignore the column corresponding to quality in most sections until the end, when we try to interpret the findings. Our first step is to read in the two tables and combine them into one. wine_red =pd.read_csv("data/wine+quality/winequality-red.csv", delimiter=";" ) wine_red['type']="red" wine_white =pd.read_csv("data/wine+quality/winequality-white.csv", delimiter=";") wine_white['type']="white" # remove spaces in column names: col_names =['fixed_acidity','volatile_acidity','citric_acid','residual_sugar', 'chlorides','free_sulfur_dioxide','total_sulfur_dioxide', 'density','pH','sulphates','alcohol','quality','type'] wine2 =pd.concat([wine_red, wine_white], ignore_index=True) wine2.columns =col_names At the end of the code chunk above, the DataFrame wine2 contains both the red and white wine data. The column names have also been edited to remove whitespaces. Here is a brief overview of the data. print(wine2.head()) fixed_acidity volatile_acidity citric_acid residual_sugar chlorides \ 0 7.4 0.70 0.00 1.9 0.076 1 7.8 0.88 0.00 2.6 0.098 2 7.8 0.76 0.04 2.3 0.092 3 11.2 0.28 0.56 1.9 0.075 54 4 7.4 0.70 0.00 1.9 0.076 free_sulfur_dioxide total_sulfur_dioxide density pH sulphates \ 0 11.0 34.0 0.9978 3.51 0.56 1 25.0 67.0 0.9968 3.20 0.68 2 15.0 54.0 0.9970 3.26 0.65 3 17.0 60.0 0.9980 3.16 0.58 4 11.0 34.0 0.9978 3.51 0.56 alcohol quality type 0 9.4 5 red 1 9.8 5 red 2 9.8 5 red 3 9.8 6 red 4 9.4 5 red 3.2 Principal Components Analysis A Principal Components Analysis (PCA) explains the covariance matrix of a set of variables through a few linear combinations of these variables. The general objectives are 1. data reduction into features that are uncorrelated with one another, 2. interpretation, and 3. visualisation. 3.2.1 Formal Set-up Suppose that we have 𝑁observations of a random vector of length 𝑝. We can represent these values in a matrix with 𝑁rows and 𝑝columns: X𝑁×𝑝 =⎡ ⎢ ⎣𝑥1,1 𝑥1,2 … 𝑥1,𝑝 ⋯ ⋯ ⋯ ⋯ 𝑥𝑁,1 𝑥𝑁,2 … 𝑥𝑁,𝑝⎤ ⎥ ⎦ Let x𝑗=[𝑥1,𝑗 𝑥2,𝑗 ⋯ 𝑥𝑁,𝑗]𝑇correspond to column 𝑗in X, for 𝑗=1,…,𝑝. We represent the mean of column 𝑗with 𝑥𝑗=1 𝑁𝑁 ∑ 𝑖=1𝑥𝑖,𝑗 The first step in a PCA is to compute the covariance matrix of the data: 𝑆𝑝×𝑝 =⎡ ⎢ ⎣𝑠2 1𝑠2 1,2 … 𝑠2 1,𝑝 ⋯ ⋯ ⋯ ⋯ 𝑠2 𝑝,1 𝑠2 𝑝,2 … 𝑠2 𝑝⎤ ⎥ ⎦ 55 where 𝑠2 𝑚,𝑛 =1 𝑁−1∑𝑁 𝑘=1(𝑥𝑘,𝑚− 𝑥𝑚)(𝑥𝑘,𝑛− 𝑥𝑛)is the sample covariance between columns 𝑚 and 𝑛of matrix 𝑋, where 1<𝑚,𝑛<𝑝. A PCA analysis yields coefficients 𝑎𝑖such that: y𝑖=𝑎𝑖,1x1+𝑎𝑖,2x2+⋯+𝑎𝑖,𝑝x𝑝, 𝑖=1,…,𝑝 In other words, y𝑖is a column vector of length 𝑁, formed from a linear combinations of the columns in the original Xmatrix. Each y𝑖is what we refer to as a principal component. From a symmetric 𝑝×𝑝matrix, we can always compute 𝑝principal components, and these vectors will be uncorrelated with each other. ĹNote However this does not help us! We have not achieved any reduction!? The value of PCA comes from the possibility that the first few principal components usually explain most of the variability in the data (∑𝑝 𝑖𝑠2 𝑖). The last few principal components typically explain little of the variability in the data. Indeed, there is some loss of information when we drop them, but the benefit is that we can (hopefully) focus on much fewer dimensions than the original 𝑝(which could be in the hundreds, even). Moreover, as these components will be uncorrelated by design, they can be used as features to solve any issues of multicollinearity in our data. Example 3.2 (PCA on wine dataset).While it is possible to extract principal components using either the covariance matrix or the correlation matrix, using the latter avoids situations where the primary principal component is simply driven by the scale of one or more columns in the original dataset. Here, we scale the first 11 columns (exclude quality and type) so that each column has mean 0 and variance 1. X_raw =wine2.iloc[:, :-2] scaler =preprocessing.StandardScaler().fit(X_raw) X_scaled =scaler.transform(X_raw) At this point, the scaler object contains all the information needed to scale a dataset, using the column means and standard deviations computed from X_raw. It is theoretically possible to extract 11 components from this 𝑋matrix. Let us proceed with that, and assess how many we should keep using a scree plot. pca_full =decomposition.PCA(n_components=11) pca_full.fit(X_scaled) n_components 11 copy True whiten False svd_solver 'auto' tol 0.0 iterated_power 'auto' n_oversamples 10 power_iteration_normalizer 'auto' 56 random_state None A scree plots the variance explained by each subsequent principal component (on the 𝑦-axis) versus the order of the principal component. It indicates how much more value there is in including the subsequent principal component. Generally, we look for an “elbow” shape to inform us of how many to keep. From below, we would probably want to retain 4 or 5 principal components. PC_values =np.arange(pca_full.n_components_) +1 plt.figure(figsize=(6,3)) plt.plot(PC_values, pca_full.explained_variance_ratio_, 'o-', linewidth=2, color='blue') plt.title('Scree Plot') plt.xlabel('Principal Component') plt.ylabel('Variance Explained'); 2 4 6 8 10 Principal Component 0.00 0.05 0.10 0.15 0.20 0.25 Variance Explained Scree Plot Figure 3.1: Scree plot To find the amount of total variance explained, we can use the following command: pca_full.explained_variance_ratio_.cumsum() array([0.2754426 , 0.50215406, 0.64364015, 0.73187216, 0.79731533, 0.85252548, 0.90008537, 0.94567722, 0.97631577, 0.99701538, 1. ]) With scikit-learn, estimated parameters will be accessible through attributes that end in an underscore. Above, we can see that the first principal component explains 27.5% of the total variation. The first two componetns explain 50.2% of the variation, and so on. 57 Example 3.4 (Wine clustering quality).Here are the silhouette scores from the two clusterings of the wine data. out =hierarchy.cut_tree(hc1, n_clusters=3).ravel() X_transformed_df['groups']=out clust.compute_silhouette_scores(hc1, X_transformed_df.iloc[:, :-2], [2,3,4]) Computing groupings for k=2 Computing score for k=2 Computing groupings for k=3 Computing score for k=3 Computing groupings for k=4 Computing score for k=4 [0.2988729940347387, 0.2549456216982205, 0.2446142652319531] The silhouette coefficient values we are obtaining are not very good. However, out of the possible values we tried, 𝐾=2seems to be the best. out =hierarchy.cut_tree(hc1, n_clusters=2).ravel() X_transformed_df['groups']=out sns.relplot(data=X_transformed_df, x='PC1', y='PC2', col='groups', hue='type', marker='o', alpha=0.3, height=3, aspect=1.2); 5 0 5 PC1 5 0 5 10 PC2 groups = 0 5 0 5 PC1 groups = 1 type red white As we can see, the groupings closely mirror the type of wine (red or white). wine2.type.groupby(X_transformed_df.groups).describe() count unique top freq groups 0 1684 2 red 1549 64 count unique top freq groups 1 4813 2 white 4763 However, notice that there are a number of white wines grouped as 0 (with most of the other reds). It would be interesting to study what qualities of these wines led to them being grouped with the reds. Also, consider what value the PCA brought to this problem. Go back and re-run the clustering algorithm with the scaled but untransformed data. Does the quality of clustering differ? 3.4 Outlier Detection One efficient way of performing outlier detection in high-dimensional datasets is to use random forests. The ensemble.IsolationForest object in scikit-learn “isolates” observations by: 1. Randomly selecting features, and then randomly selecting a split value between the maximum and minimum values of the selected feature to form a decision tree (see Section 8.6.1). • In each tree, the number of splittings required to isolate a sample is equivalent to the path length from the root node to the terminating node. 2. Repeating step 1 to create a forest of trees (see Section 8.6.3). For each observation, the average path length, over the forest of random trees, is a measure of how anomalous it is. Random partitioning produces noticeably shorter paths for anomalies. Hence, when a forest of random trees collectively produce shorter path lengths for particular samples, they are highly likely to be anomalies. Example 3.5 (Isolation forest with taiwan dataset).Let us apply this technique to the Taiwan real estate dataset from the regression topic. We shall see much more of this dataset in Section 5.1. re2 =pd.read_csv("data/taiwan_dataset.csv") X_re =re2.loc[:, ['trans_date','house_age','dist_MRT','num_stores', 'Xs','Ys','price']] X_re_scaled =preprocessing.StandardScaler().fit_transform(X_re) After reading in and scaling the data, we fit the IsolationForest estimator. clf =IsolationForest(max_samples=300, max_features=2, contamination=0.01, random_state=503) clf.fit(X_re_scaled) n_estimators 100 max_samples 300 contamination 0.01 max_features 2 65 bootstrap False n_jobs None random_state 503 verbose 0 warm_start False Here is a brief explanation of the arguments in the call to IsolationForest: •contamination factor is the proportion of outliers that we expect to see in the dataset. •max_features corresponds to the number of features to be drawn for each base estimator •max_samples is the number of samples (observations) to draw from the original dataset for each base estimator. id_outliers =pd.Series(clf.predict(X_re_scaled)) id_outliers.value_counts() 1 409 -1 5 Name: count, dtype: int64 There are 5 points that have been identified as outliers (coded as -1). As analysts, we should do our best to understand what property, or combination of features, led to this. re2['outliers']=id_outliers outliers -1 1 house_age 25.900 17.612 dist_MRT 6042.907 1023.262 num_stores 1.000 4.132 price 14.920 38.262 Xs -5.424 0.066 Ys -1.982 0.024 We can see that the outliers are very different from the remaining points: the distance to MRT, X-coordinates, and price are all very different. Having identified these points, our job as analysts is to interpret the differences. ĹNote If we drop the location parameters, would different points be identified as outliers? 66 3.5 Visualisation 3.5.1 MDS Multidimensional Scaling is a technique for visualising high-dimensional data. Just like in hierarchical clustering, we begin with a square matrix consisting of all pairwise dissimilarities 𝑑(𝑥𝑖,𝑥𝑗)between our 𝑁high-dimensional vectors. With a choice of 𝑘, we seek values 𝑧1,𝑧2,…,𝑧𝑁∈ℝ𝑘such that the following function is minimised: 𝑆(𝑧1,…,𝑧𝑁)=[∑ 𝑖≠𝑗(𝑑(𝑥𝑖,𝑥𝑗)−||𝑧𝑖−𝑧𝑗||)2]1/2 Thus the goal of MDS is to find a lower dimensional set of vectors whose pairwise Euclidean distances are as close as possible to the dissimilarity matrix of the original vectors. MDS is not the same as Principal Component Analysis (PCA): • PCA maximises variance, orthogonal to earlier components. • Principal components are ordered; MDS are not. • Principal components are linear combinations of the original vectors; MDS output is not. Example 3.6 (MDS on disease symptoms).The dataset disease.csv contains a list of symptoms that were reported for a set diseases. Each row in the dataframe corresponds to a particular disease, while each binary column indicates whether that particular symptom was frequently present for this disease. disease =pd.read_csv("data/disease.csv") Our goal is to visualise the 41 diseases - diseases with “similar” symptoms should be plotted “close” to one another. However, the data begets the natural question: How can we define dissimilarity between the symptom lists of two diseases? For this purpose, we shall use the Jaccard similarity index. For two disease symptom sets 𝐴 and 𝐵, the Jaccard dissimilarity is defined to be 𝐽=1−|𝐴∩𝐵| |𝐴∪𝐵| If there are no common symptoms between the two diseases, then the intersection between the two sets would be the empty set. In that situation, 𝐽would take on the maximum value of 1. If the two symptom lists are identical, then 𝐽takes on the smallest possible value of 0. In this next chunk of code, we concatenate all the symptoms into a single string, separated by commas. 67 disease_names =disease.disease.to_list() symptoms =disease.columns.to_list()[:-1] X=disease.iloc[:, 0:-1].to_numpy() symptom_text =[] for iin range(0, X.shape[0]): symptom_text.append(','.join([symptoms[x] for xin np.where(X[i] == 1)[0]])) disease['symptom_text']=symptom_text For instance, here are the symptoms for GERD and heart attack. disease.loc[disease.disease.isin(['GERD', 'Heart attack']), 'symptom_text'].to_list() ['stomach_pain,acidity,ulcers_on_tongue,vomiting,cough,chest_pain', 'vomiting,chest_pain,breathlessness,sweating'] Based on the set of symptoms for GERD and Heart attack, we have the following Jaccard calculation: 𝐽=1−2 8=0.75 Now we turn to the MDS transformation. embedding =MDS(n_components=2, normalized_stress='auto', dissimilarity='precomputed', n_init=4, random_state=42, max_iter=500, verbose=0) # pdist2 is 41x41 pdist2 =pairwise_distances(X!=0, metric='jaccard') X_transformed =embedding.fit_transform(pdist2) X_transformed_df =pd.DataFrame(X_transformed, columns=['X','Y']) X_transformed_df['disease']=disease_names In this pdf version, we display a static image. However, the online version contains an interactive plot. plt.figure(figsize=(8,4)) ax =sns.scatterplot( data=X_transformed_df, x="X", y="Y" ) for _, r in X_transformed_df.iterrows(): ax.text(r["X"], r["Y"], r["disease"], ha="center", va="bottom", fontsize=8) 68 0.6 0.4 0.2 0.0 0.2 0.4 0.6 0.8 X 0.8 0.6 0.4 0.2 0.0 0.2 0.4 0.6 0.8 Y Fungal infection Allergy GERD Chronic cholestasis Drug reaction Peptic ulcer disease AIDS Diabetes Gastroenteritis Bronchial asthma Hypertension Migraine Cervical spondylosis Paralysis (brain hemorrhage) Jaundice Malaria Chicken pox Dengue Typhoid Hepatitis A Hepatitis B Hepatitis C Hepatitis D Hepatitis E Alcoholic hepatitis Tuberculosis Common Cold Pneumonia Dimorphic hemmorhoids(piles) Heart attack Varicose veins Hypothyroidism Hyperthyroidism Hypoglycemia Osteoarthristis Arthritis Vertigo Acne Urinary tract infection Psoriasis Impetigo To verify if the plot makes intuitive sense, we should inspect the points that occur nearby to one another. For instance, since we observe Typhoid and Malaria to be close on the plot, we can retrieve their actual symptoms with: for xin disease.loc[disease.disease.isin(['Typhoid','Malaria']), 'symptom_text'].to_list(): pprint(x) 'chills,vomiting,nausea,high_fever,diarrhoea,headache,sweating,muscle_pain' 'chills,vomiting,nausea,abdominal_pain,high_fever,fatigue,diarrhoea,headache,constipation,toxic_look_(typhos),belly_pain' It is important to remember that above, we requested the MDS algorithm to return us 𝑧coordinates in R2two-dimensional Euclidean space, that replicate the dissimilarity matrix from the higher-dimensional data. This was for convenience, since the two dimensional plane is easier to plot. One problem is that, due to information loss, it is possible that in 2D, two points appear near, but in fact, they may be far apart on a third dimension. With plotly, we can make interactive 3d plots, which can go some way to alleviating this problem. ÁWarning Stay away from non-interactive 3d plots! For an interaactive version of the 3D plot, please visit the online version of the textbook. 3.5.2 t-SNE t-SNE is an fast, iterative algorithm for visualising high-dimensional data. Once again, suppose we have 𝑁data points 𝑥𝑖∈R𝑝. We would like to choose 𝑁map points 𝑦𝑖∈R2to represent them. Here is how the algorithm works: 69 1. Compute pairwise similarity between the data points, using a Gaussian (Normal) kernel. 2. Iteratively update map points so that their pairwise similarity is as close as possible to the original data points. The innovation of these algorithm is that the similarity between map points is computed using a𝑡-distribution instead of Gaussian. The 𝑡-distribution has fatter tails than the Normal. This ensures that data points that are not close in R𝐷are pushed apart in the map points space. There are a couple of important parameters in this algorithm. 1. Perplexity: This is a parameter that provides a guide on how many neighbours a point has. It is recommended to try different values between 5 and 50 and assess if results are meaningful and consistent. 2. Number of iterations: This is the number of adjustments to make to the map points before stopping. The difference between map point similarity and data point similarity is measured using Kullback-Leibler Divergence. When this no longer drops quickly, we can stop the t-SNE algorithm. Although you will find t-SNE used in many analyses, you should be aware of certain caveats when using it. Please take a look at the link in the references for plots that expound on the following points. 1. t-SNE is for visualisation and exploration, not for clustering. 2. Distances in the map space are not reflective of true distances between datapoints. 3. t-SNE preserves small pairwise distances so that local relationships are preserved. 4. In the map space, distances between clusters (between far-away points) might not mean much (unlike MDS). 5. Make sure you try with different perplexity values, and check that the algorithm has converged. Example 3.7 (Twitter dataset).The UCI Machine Learning repository contains tweets pertaining to health news from more than 15 major news agencies in 2015. In this example, we are going to encode each BBC tweet in to a numerical vector of length 384. Then we are going to visualise it using t-SNE. For more details on how words are converted into vectors, take a look at Section 4.5. In the next code chunk, we use a neural language model to convert each tweet into the vector of length 384. The model is one of several models provided and maintained by Hugging Face model =SentenceTransformer('sentence-transformers/all-MiniLM-L12-v2') #sentences = ["This is an example sentence. can we move on?", # "Each sentence is converted"] #embeddings = model.encode(sentences) #embeddings.shape Now we turn to the dataset of tweets. bbchealth_df =pd.read_table('data/health+news+in+twitter/Health-Tweets/bbchealth.txt', delimiter='|', names=['id','datetime','tweet']) bbchealth_tweets =bbchealth_df.tweet print(bbchealth_tweets[10]) 70 Have GP services got worse? http://bbc.in/1Ci5c22 Above, we have an example of a tweet. We are going to strip off the URL at the end before encoding each (short) sentence. t1 =bbchealth_tweets.str.replace(' http.*$','', regex=True) t2 =t1.str.replace('^VIDEO:','', regex=True) t2_l =t2.to_list() embeddings =model.encode(t2_l) Example 3.8 (Twitter dataset t-SNE output).The following code generates an interactive plot based on the t-SNE visualisation. tsne1 =TSNE(n_components=2, init="random", perplexity=10, verbose=0, random_state=43, max_iter=5000) X_transformed2 =tsne1.fit_transform(embeddings) df2 =pd.DataFrame(X_transformed2, columns=['x','y']) df2['labels']=t2_l Please refer to the online version of the text for an interactive plot. plt.figure(figsize=(8,4)) sns.scatterplot( data=df2, x="x", y="y") 150 100 50 0 50 100 x 150 100 50 0 50 100 150 y For another application of t-SNE visualisation, refer to Section 4.6. 71 3.6 References 3.6.1 Website references 1. Clustering performance evaluation There are many other ways of assessing what the optimal number of clusters should be. 2. More information on isolation forests 3. Wikipedia entry on Jaccard Similarity (or index): This page introduces variants of the index, which may be relevant when one has counts of the number of times each item appears in the set. 4. Caveats when using t-SNE 5. DataCamp on t-SNE: This tutorial consists of a worked-through example on a churn dataset. It includes a comparison with PCA. 3.6.2 Video references 1. t-SNE, explained by Josh Starmer 2. A general video on high-dimensional space: A nice explainer from Google. 72 4 Natural Language Processing 4.1 Introduction In our world, Natural Language Processing (NLP) is used in several scenarios. For example, • mobile phones and personal computers support predictive text. • web search engines give access to information locked up in unstructured text; • machine translation allows us to understand texts written in languages that we do not know; • text analysis enables us to detect sentiment in tweets and blogs. But as we begin to explore NLP, we realise that it is an extremely difficult subject. Here are some examples to note: 1. Some words mean different things in different contexts, but as humans, we know which meaning is being used. • He served the dish. 2. In the following two sentences, the word “by” has different meanings: • The lost children were found by the lake. • The lost children were found by the search party. 3. In the following cases, we (humans) can resolve what “they” is referring to, but it is not easy to generate a simple rule that a computer can follow. • The thieves stole the paintings. They were subsequently recovered. • The thieves stole the paintings. They were subsequently arrested. 4. How can we get a computer to understand that the following tweet carries a negative sentiment? • “Wow. Great job st@rbuck’s. Best cup of coffee ever.” ĹNote Can you catch all three jokes in the movie clip below? � https://youtu.be/NfN_gcjGoJo import numpy as np import pandas as pd from itables import show import pprint import gensim 73 vectorizer2 =TfidfVectorizer(stop_words='english', norm=None) X2 =vectorizer2.fit_transform(raw_docs) print(pd.DataFrame(X2.A, columns=list(vectorizer2.get_feature_names_out())).iloc[:, :10].round(3)) afraid basic comes counting data documents examples interesting \ 0 0.000 1.288 0.000 0.000 0.000 0.000 0.000 0.000 1 1.288 0.000 0.000 0.000 0.000 0.000 0.000 1.693 2 1.288 2.575 1.693 1.693 1.693 1.693 1.693 0.000 just larger 0 0.000 0.000 1 0.000 0.000 2 1.693 1.693 The final step normalises the weights within each document. vectorizer3 =TfidfVectorizer(stop_words='english') X3 =vectorizer3.fit_transform(raw_docs) print(pd.DataFrame(X3.A, columns=list(vectorizer3.get_feature_names_out())).iloc[:, :10].round(3)) afraid basic comes counting data documents examples interesting \ 0 0.000 0.577 0.000 0.000 0.000 0.000 0.000 0.000 1 0.474 0.000 0.000 0.000 0.000 0.000 0.000 0.623 2 0.170 0.340 0.223 0.223 0.223 0.223 0.223 0.000 just larger 0 0.000 0.000 1 0.000 0.000 2 0.223 0.223 The above matrix is known as a document-term matrix, since the columns are defined by terms, and each row is a document. At this point, we can use each row as a vector representation of each document. If necessary, for this corpus, we could even represent each term using its corresponding column. Note that some books/software use a slightly different convention - they may work with the term-document matrix. However, the idea is the same. Take a look at the following termdocument matrix, assembled from the complete works of Shakespeare: Figure 4.1: Jurafsky and Martin (2025) 80 If we intend to represent each document as a numeric vector, the columns, highlighted by the red boxes, would be a natural choice. Suppose we only focus on the coordinates corresponding to the words battle and fool. Then a visualisation of the documents would look like this: Figure 4.2: Jurafsky and Martin (2025) Visually, it is easy to tell that “Henry V” and “Julius Caesar” are similar (they point in the same direction) as opposed to “As You Like It” and “Twelfth Night”. But it is also easy to see why - the former two contain similar high counts of battle compared to the latter two, which are comedies. Tf-idf are a normalised version of the above raw counts; they provide a numerical representation of documents, adjusting for document length and words that are common across all documents in a corpus. 4.5.2 Cosine similarity In order to quantify the similarity (or nearness) of vector representations in NLP, the common method used is cosine similarity. Suppose that we have a vector representation of two documents vand w. If the vocabulary size is 𝑁, then each of the vectors is of length 𝑁. Since we are dealing with counts the coordinate values of each vector will be non-negative. We use the angle 𝜃between the vectors as a measure of their similarity: cos 𝜃= ∑𝑁 𝑖=1𝑣𝑖𝑤𝑖 √∑𝑁 𝑖=1𝑣2 𝑖√∑𝑁 𝑖=1𝑤2 𝑖 Geometrically, cosine similarity measures the size of the angle between vectors: 81 Figure 4.3: Jurafsky and Martin (2025) 4.5.3 Dense Embeddings One of the drawbacks of sparse vectors is that they are very long (the length of the vocabulary), and most entries in the vector will be 0. As a result, researchers worked on methods that would pack the information in the vectors into shorter ones. Instead of working on representations of the documents, the methods aimed to create representations of each token (or word) in the vocabulary. These are referred to as embeddings. Here, we shall discuss word2vec (Mikolov et al. (2013)), but take note that there are others. GLoVe (Pennington, Socher, and Manning (2014)) was invented soon after, but the most common embeddings used today arise from Deep Learning models. The most widely used version is BERT (see the video references below, as well as Devlin et al. (2019)). The approach in word2vec deviates considerably from tf-idf, in that the goal is to obtain a numeric representation of a word, in the context of it’s surrounding words. Consider this statement: 13% of the United States population eats pizza on any given day. Mozzarella is commonly used on pizza, with the highest quality mozzarella from Naples. In Italy, pizza served in formal > settings is eaten with a fork and knife. The words eats,served and mozzarella appear close to pizza. Hence another word that appears in similar contexts, should be similar to pizza. Examples could be certain baked dishes or even salad. To achieve such a representation, word2vec runs a self-supervised algorithm, with two tasks: 1. Primary task: To “learn” a numeric vector that represents each word. 2. Pretext task (stepping stone): To train a classifier that, when given a word 𝑤, predicts nearby context words 𝑐. Self-supervised algorithms differ from supervised algorithms in that there are no labels that need to be created. The pre-text task trains a model to perform predictions, based on a sliding window context: Starting with an initial random vector for each word, the algorithm updates the vectors as it proceeds through the corpus, finally ending up with an embedding for each word that reflects its semantic value, based on neighbouring words. In NLP, the quality of an embedding can be evaluated using an analogy task: 82 Given X, Y and Z, find W such that W is related to to Z in the same way that X is related to Y. For instance, if we are given the pair man:king, and the word woman, then the embedding should return queen, since woman:queen in the same way that man is related to king. Geometrically, the answer to the analogy is obtained by adding (king - man) to woman. The nearest embedding to the result, is returned as the answer. On the left are examples of the types of analogy pairs that word2vec is able to solve, while on the right, we have visualisations of GLoVe. (a) word2vec (a) GloVE Example 4.5 (Glove dense embeddings).Here’s how we can use gensim code to conduct the analogy task. First, we load 400,000 GLoVe vectors, each representing a different word. Each vector is of length 100. word_vectors =api.load("glove-wiki-gigaword-100") The word_vectors object is similar to a dictionary. For instance, word_vectors['woman'] will return the vector corresponding to “woman”. The following code will run the analogy task. It returns the answer to: Man is to king as woman is to ______ . 83 # Check the "most similar words", using the default "cosine similarity" measure. result =word_vectors.most_similar(positive=['woman','king'], negative=['man']) most_similar_key, similarity =result[0]# look at the first match print(f"{most_similar_key}:{similarity:.4f}") queen: 0.7699 The object also contains methods to identify which word in a group is most dissimilar to the rest. The following code identifies “cereal” as the odd word out. print(word_vectors.doesnt_match("breakfast cereal dinner lunch".split())) cereal 4.6 Visualisation with t-SNE When compared with sparse embeddings, dense embeddings are compact. However, a more important difference is that dense vectors contain the semantic meaning of words. This means that vectors that are close to each other are similar in meaning. Let us use t-SNE to visualise the GloVe embeddings. There are a total of 400,000 vectors in the embedding. Even with t-SNE that will be difficult to make sense of. Hence for now, we simply visualise the first 100 most common words. The following code extracts the first 100 word vectors and applies the t-SNE transformation to them. Please refer to Section 3.5.2 for more details. nn =1000 X=np.zeros((nn, 100)) for ii in np.arange(nn): X[ii,] =word_vectors.get_vector(ii) labels =pd.Series([word_vectors.index_to_key[x] for xin np.arange(nn)]) tsne1 =manifold.TSNE(n_components=2, init="random", perplexity=10, metric='cosine', verbose=0, max_iter=5000, random_state=222) X_transformed2 =tsne1.fit_transform(X) You should get the same plot as us since we have set the same seed at the start of the cell, and when we initialise the transformer. Explore the resulting plot - notice how months of the year appear close together at the bottom left. Around the left as well, the calendar years appear as a group. Please refer to the online version of the text for an interactive plot, with text labels. plt.figure(figsize=(8,4)) sns.scatterplot(data=df2, x="x", y="y") 84 100 75 50 25 0 25 50 75 x 60 40 20 0 20 40 60 80 y ĹNote Would we be able to make such a plot using tf-idf? Why or why not? 4.7 Neural Language Models Language is complex. It is incredible how we can understand such long paragraphs of texts with such ease. We somehow seem to have learnt complicated sets of grammar and syntax just by listening to others speak. To get a machine to learn language has not been easy. It is only recently that Large Language Models such as chatGPT have demonstrated that it is possible for machines to converse with humans just as we do to one another. Neural Models (or deep learning models) have been the key to this. In this subsection, we provide a very brief overview of their characteristics that allow them to achieve impressive performance on a range of language-related tasks. The basic unit of a neural model is the neural unit (on the left). It consists of weights and a non-linear activation function. Given an input vector, the weights are multiplied by the input vector, summed and then fed through the activation function to generate an output. Neural models are made up of many neural units, organised into layers. The first neural models were Feed-Forward Networks. Due to the virtue of being able to incorporate many parameters, and due to semi-supervised learning, they were already a huge improvement over earlier models. Here is a simple set up, with one hidden layer for training a language model (used to predict the next word). It can also be used to learn embeddings. 85 (a) Jurafsky and Martin (2025) (a) Neuron figure from https://en.wikipedia.org/w iki/Neuron Figure 4.8: Jurafsky and Martin (2025), FFN The next evolution in neural models was the ability to incorporate words in the recent history. For humans, this comes naturally. For instance, we know that this is grammatically correct: The flights the airline was cancelling were full. For neural models to have this ability, it was necessary to incorporate the hidden layers from recent words when processing the current word. Recurrent Neural Networks (RNNs) and LongShort Term Memory (LSTM) networks had these features, but they were very slow to train. The major breakthrough came with the invention of the transformer architecture. The self-attention layer of these networks gave a word access to all preceding words in the training window, instead of just one. Most importantly. the training of these networks could be parallelised! 86 Figure 4.9: Jurafsky and Martin (2025) Here are some more examples where transformers excel: The keys to the cabinet are on the table. The chicken crossed the road because it wanted to get to the other side. I walked along the pond, and noticed that one of the trees along the bank had fallen into the water after the storm. In the final sentence, the word bank has two meanings - how will a model know to decide the correct one? With transformers, because the full context of a word is captured along with it, it is possible to perform this disambiguation. Figure 4.10: Jurafsky and Martin (2025) 4.8 Applications Hugging Face has spent a considerable effort to make Neural Language Models accessible and available to all with minimal coding. For starters, they have ensured that all their models are described in a standardised manner with model cards. Here is an example of a model card for BERT. 87 Moreover, they have developed easy to use pipelines. For NLP, the following tasks have mature pipelines: • feature-extraction (obtaining the embedding of a text) • ner • question-answering • sentiment-analysis • summarization • text-generation • translation, and • zero-shot-classification. 4.8.1 Sentiment Analysis In this subsection, we shall utilise one of their sentiment analysis models on the wine reviews dataset. This is a transformer-based neural language model (BERT) that has been fine-tuned with data labelled with sentiments. All we have to do is feed in the sentence, and we will obtain a confidence score, and a sentiment label. classifier =pipeline("sentiment-analysis", model="distilbert/distilbert-base-uncased-finetuned-sst-2-english") classifier(["I love this course!","I absolutely detest this course."]) Device set to use cpu [{'label': 'POSITIVE', 'score': 0.9998835325241089}, {'label': 'NEGATIVE', 'score': 0.9973570704460144}] Example 4.6 (Wine reviews sentiments).The number of reviews we have is close to 120,000. Hence, computing the sentiments for each and every one will take a long time. Instead, we shall compute the sentiments for a sample (of size 20, where possible) from each variety of wine. The following snippet samples 20 reviews from each wine type. tmp_df =wine_reviews.head(0).copy() for x,vv in wine_reviews.groupby(wine_reviews.variety): grp_len =vv.shape[0] if(grp_len >= 20): vv =vv.sample(n=20, random_state=99) tmp_df =pd.concat([tmp_df, vv], ignore_index=True) review_list =list(tmp_df.description) The next snippet computes the sentiment scores for those sampled reviews. 88 tmp_df['score']=0.00 tmp_df['label']='' for i,rr in enumerate(review_list): tmp =classifier(rr)[0] tmp_df.loc[i, 'score']=tmp['score'] tmp_df.loc[i, 'label']=tmp['label'] This next snippet tabulates the sentiment classifications for the wine types. sent_counts =pd.crosstab(tmp_df.variety, tmp_df.label, margins=True) sent_counts['proportion']=sent_counts.POSITIVE/sent_counts.All sent_counts.head() label NEGATIVE POSITIVE All proportion variety Abouriou 0 3 3 1.000000 Agiorgitiko 0 20 20 1.000000 Aglianico 0 20 20 1.000000 Aidani 0 1 1 1.000000 Airen 1 2 3 0.666667 These are the reviews for one of the varieties that had a proportion of positive reviews close to 50%. for xin wine_reviews[wine_reviews.variety == 'Tempranillo Blanco'].description.values: pp.pprint(x) ("Gold in color and lightly oxidized on the nose, and it's still young. Smells " 'heavy and creamy, like hay. Feels flat, with pickled flavors and mealy apple ' 'on the finish. Runs plump, sweet and seems like an imposter for Chardonnay.') ('Oily, stalky, bready aromas are a bit tired. This has a chunky feel offset ' 'by citric acidity. Briny, salty flavors of citrus fruits and lees are ' "lasting. For varietal Tempranillo Blanco, this isn't bad.") ('Maderized in color, this wine has a yeasty, creamy nose with baked ' "white-fruit aromas and caramel. It's OK in feel, with pickled, mildly briny " 'flavors of apple and apricot. The finish is showing some oxidization, ' 'leading to a chunky, fleshy feel.') ("Waxy peach aromas seem slightly oxidized. It's round and citrusy on the " 'palate, but in a monotone way that fades to pithy white fruits and mealy ' "citrus. Shows some flashes of uniqueness and class; mostly it's wayward and " 'slightly bitter.') ('Forget the high price on this Tempranillo Blanco. Looking at the wine alone, ' "it's briny and stalky on the nose, with wiry lemon-like acids that push sour " "orange flavors. Overall it's monotone, briny and citrusy.") ("A maderized color is apropos for the wine's fully mature, nutty nose. This " 'is big and cidery feeling, with apple and orange flavors. A finish of ' 89 Figure 4.13: Sun et al. (2021) Figure 4.14: Clark et al. (2019) 96 neural models (11:37) 2. Introduction to RNN 3. Introduction to BERT 4.10.2 Website References 1. Hugging Face course on transformers 2. Gensim documentation: Contains tutorials as well. 3. Using sklearn to perform LDA: We can also use scikit-learn to perform LDA. 4. Visualising LDA: Contains sample notebooks for the visualisation. 97 5 Linear Regression 5.1 Introduction Regression analysis is a technique for investigating and modeling the relationship between variables like X and Y. Here are some examples: 1. Within a country, we may wish to use per capita income (X) to estimate the life expectancy (Y) of residents. 2. We may wish to use the size of a crab claw (X) to estimate the closing force that it can exert (Y). 3. We may wish to use the height of a person (X) to estimate their weight (Y). In all the above cases, we refer to 𝑋as the explanatory or independent variable. It is also sometimes referred to as a predictor.𝑌is referred to as the response or dependent variable. In this topic, we shall first introduce the case of simple linear regression, where we model the 𝑌 on a single 𝑋. In later sections, we shall model the 𝑌on multiple 𝑋’s. This latter technique is referred to as multiple linear regression. Regression models are used for two primary purposes: 1. To understand how certain explanatory variables affect the response variable. This aim is typically known as estimation, since the primary focus is on estimating the unknown parameters of the model. 2. To predict the response variable for new values of the explanatory variables. This is referred to as prediction. In this topic, we shall focus on the estimation aim, since prediction models require a paradigm of their own, and are best learnt alongside a larger suite of models e.g. decision trees, support vector machines, etc. We shall cover prediction in the topic of supervised learning Section 8.1. import pandas as pd import numpy as np import matplotlib.pyplot as plt import seaborn as sns import folium import geopandas from itables import show import statsmodels.api as sm from statsmodels.formula.api import ols from scipy import stats 98 Example 5.1 (Taiwan real estate).For this tutorial, we shall work with a data set from the UCI machine learning repository. It contains real estate prices in the Xindian district of Taiwan. Our goal is to answer the following question: How well can we explain real-estate prices in Taiwan? Here is a brief description of the columns in the dataset: •trans_date: The date of the transaction. As you can see, this has been coded to be a numerical value, so 2013.5 refers to June 2013. •house_age: Age of the house in years. •dist_MRT: Distance the the nearest MRT (in metres) •num_stores: Number of convenience stores within walking distance •lat,long: Latitude and longitude •price: House price per unit area (10,000 New Taiwan Dollars per Ping, which is about 3.3𝑚2) •X,Y,Xs,Ys: Projected coordinates re2 =pd.read_csv("data/taiwan_dataset.csv") print(re2.head()) id trans_date house_age dist_MRT num_stores lat long \ 0 1 2012.916667 32.0 84.87882 10 24.98298 121.54024 1 2 2012.916667 19.5 306.59470 9 24.98034 121.53951 2 3 2013.583333 13.3 561.98450 5 24.98746 121.54391 3 4 2013.500000 13.3 561.98450 5 24.98746 121.54391 4 5 2012.833333 5.0 390.56840 5 24.97937 121.54245 price X Y Xs Ys 0 37.9 506501.554580 2.766295e+06 0.666075 1.559003 1 42.2 506433.292007 2.766001e+06 0.597812 1.265025 2 47.3 506862.973688 2.766798e+06 1.027494 2.062484 3 54.8 506862.973688 2.766798e+06 1.027494 2.062484 4 43.1 506732.307872 2.765899e+06 0.896828 1.163084 Here is a scatter plot of the points in the dataset. ax =sns.scatterplot(data=re2, x='X', y='Y', hue='price', size='price') ax.set_title("Taiwan real estate data"); 99 500000 502000 504000 506000 508000 X 2.762 2.764 2.766 2.768 2.770 Y 1e6 Taiwan real estate data price 20 40 60 80 100 Let us first explore the dataset with Python. How do this plots aid in your understanding of the dataset? plt.figure(figsize=(10,4)); ax=plt.subplot(131) re2.plot(x='dist_MRT', y='price', kind='scatter', ax=ax, title='distance to MRT') ax=plt.subplot(132) re2.plot(x='house_age', y='price', kind='scatter', ax=ax, title='House age') ax=plt.subplot(133) z=re2.num_stores.unique() z.sort() tmp2 =np.array([re2.price[re2.num_stores == x].to_numpy() for xin z], dtype=object) ax.boxplot(tmp2, tick_labels=z);ax.set_title('Num. of stores'); 100 0 2000 4000 6000 dist_MRT 20 40 60 80 100 120 price distance to MRT 0 10 20 30 40 house_age 20 40 60 80 100 120 price House age 012345678910 20 40 60 80 100 120 Num. of stores 5.2 Simple Linear Regression 5.2.1 Formal Set-up The simple linear regression model is applicable when we have observations (𝑋𝑖,𝑌𝑖)for 𝑛 individuals. For now, let’s assume both the 𝑋and 𝑌variables are quantitative. The simple linear regression model is given by 𝑌𝑖=𝛽0+𝛽1𝑋𝑖+𝑒𝑖 where •𝛽0is intercept term, •𝛽1is the slope, and •𝑒𝑖is an error term, specific to each individual in the dataset. 𝛽0and 𝛽1are unknown constants that need to be estimated from the data. There is an implicit assumption in the formulation of the model that there is a linear relationship between 𝑌𝑖and 𝑋𝑖. In terms of distributions, we assume that the 𝑒𝑖are i.i.d Normal. 𝑒𝑖∼𝑁(0,𝜎2),𝑖=1…,𝑛 The constant variance assumption is also referred to as homoscedascity (homo-skee-das-city). The validity of the above assumptions will have to be checked after the model is fitted. All in all, the assumptions imply that: 1. 𝐸(𝑌𝑖|𝑋𝑖)=𝛽0+𝛽1𝑋𝑖, for 𝑖=1,…,𝑛. 2. 𝑉𝑎𝑟(𝑌𝑖|𝑋𝑖)=𝑉𝑎𝑟(𝑒𝑖)=𝜎2, for 𝑖=1,…,𝑛. 3. The 𝑌𝑖are independent. 4. The 𝑌𝑖’s are Normally distributed. 101 5.2.2 Estimation Before deploying or using the model, we need to estimate optimal values to use for the unknown 𝛽0and 𝛽1. We shall introduce the method of Ordinary Least Squares (OLS) for the estimation. Let us define the error Sum of Squares to be 𝑆𝑆𝐸=𝑆(𝛽0,𝛽1)= 𝑛 ∑ 𝑖=1(𝑌𝑖−𝛽0−𝛽1𝑋𝑖)2 Then the OLS estimates of 𝛽0and 𝛽1are given by arg min 𝛽0,𝛽1 𝑛 ∑ 𝑖=1(𝑌𝑖−𝛽0−𝛽1𝑋𝑖)2 The minimisation above can be carried out analytically, by taking partial derivative with respect to the two parameters and setting them to 0. 𝜕𝑆 𝜕𝛽0= −2 𝑛 ∑ 𝑖=1(𝑌𝑖−𝛽0−𝛽1𝑋𝑖)=0 𝜕𝑆 𝜕𝛽1= −2 𝑛 ∑ 𝑖=1𝑋𝑖(𝑌𝑖−𝛽0−𝛽1𝑋𝑖)=0 Solving and simplifying, we arrive at the following:  𝛽1=∑𝑛 𝑖=1(𝑋𝑖− 𝑋)(𝑌𝑖− 𝑌) ∑𝑛 𝑖=1(𝑋𝑖− 𝑋)2  𝛽0= 𝑌−  𝛽0 𝑋 where  𝑌 =(1/𝑛)∑𝑌𝑖and  𝑋=(1/𝑛)∑𝑋𝑖. If we define the following sums: 𝑆𝑋𝑌 =𝑛 ∑ 𝑖=1𝑋𝑖𝑌𝑖−(∑𝑛 𝑖=1𝑋𝑖)(∑𝑛 𝑖=1𝑌𝑖) 𝑛 𝑆𝑋𝑋 =𝑛 ∑ 𝑖=1𝑋2 𝑖−(∑𝑛 𝑖=1𝑋𝑖)2 𝑛 then a form convenient for computation of  𝛽1is  𝛽1=𝑆𝑋𝑌 𝑆𝑋𝑋 Once we have the estimates, we can use the estimated model to compute fitted values for each observation, corresponding to our best guess of the mean of the distributions from which the observations arose:  𝑌𝑖= 𝛽0+ 𝛽1𝑋𝑖, 𝑖=1,…,𝑛 102 As always, we can form residuals as the deviations from fitted values. 𝑟𝑖=𝑌𝑖− 𝑌𝑖 Residuals are our best guess at the unobserved error terms 𝑒𝑖. Squaring the residuals and summing over all observations, we can arrive at the following decomposition, which is very similar to the one in the ANOVA model: 𝑛 ∑ 𝑖=1(𝑌𝑖− 𝑌)2 ⏟⏟⏟⏟⏟ 𝑆𝑆𝑇 =𝑛 ∑ 𝑖=1(𝑌𝑖− 𝑌𝑖)2 ⏟⏟⏟⏟⏟ 𝑆𝑆𝑅𝑒𝑠 +𝑛 ∑ 𝑖=1( 𝑌𝑖− 𝑌)2 ⏟⏟⏟⏟⏟ 𝑆𝑆𝑅𝑒𝑔 where •𝑆𝑆𝑇is known as the total sum of squares. •𝑆𝑆𝑅𝑒𝑠 is known as the residual sum of squares. •𝑆𝑆𝑅𝑒𝑔 is known as the regression sum of squares. In our model, recall that we had assumed equal variance for all our observations. We can estimate 𝜎2with  𝜎2=𝑆𝑆𝑅𝑒𝑠 𝑛−2 Our distributional assumptions lead to the following for our estimates  𝛽0and  𝛽1:  𝛽0∼ 𝑁(𝛽0,𝜎2(1/𝑛+  𝑋2/𝑆𝑋𝑋)) (5.1)  𝛽1∼ 𝑁(𝛽1,𝜎2/𝑆𝑋𝑋)(5.2) The above are used to construct confidence intervals for 𝛽0and 𝛽1, based on 𝑡-distributions. 5.2.3 Hypothesis Test for Model Significance The first test that we introduce here is to test if the coefficient 𝛽1is significantly different from 0. It is essentially a test of whether it was worthwhile to use a simple linear regression, instead of a simple mean to represent the data. The null and alternative hypotheses are: 𝐻0∶ 𝛽1=0 𝐻1∶ 𝛽1≠0 The test statistic is 𝐹0=𝑆𝑆𝑅𝑒𝑔/1 𝑆𝑆𝑅𝑒𝑠/(𝑛−2) Under the null hypothesis, 𝐹0∼𝐹1,𝑛−2. 103 It is also possible to perform this same test as a 𝑡-test, using the result earlier. The statement of the hypotheses is equivalent to the 𝐹-test. The test statistic 𝑇0= 𝛽1 √ 𝜎2/𝑆𝑋𝑋 Under 𝐻0, the distribution of 𝑇0is 𝑡𝑛−2. This 𝑡-test and the earlier 𝐹-test in this section are identical. It can be proved that 𝐹0=𝑇2 0; the obtained 𝑝-values will be identical. 5.2.4 Coefficient of Determination, 𝑅2 The coefficient of determination 𝑅2is defined as 𝑅2=1−𝑆𝑆𝑅𝑒𝑠 𝑆𝑆𝑇=𝑆𝑆𝑅𝑒𝑔 𝑆𝑆𝑇 It can be interpreted as the proportion of variation in 𝑌𝑖, explained by the inclusion of 𝑋𝑖. Since 0≤𝑆𝑆𝑅𝑒𝑠 ≤𝑆𝑆𝑇, we can easily prove that 0≤𝑅2≤1. The larger the value of 𝑅2is, the better the model is. When we get to the case of multiple linear regression, take note that simply including more variables in the model will increase 𝑅2. This is undesirable; it is preferable to have a parsimonious model that explains the response variable well. Example 5.2 (Price vs. house age).As a first model, we fit price (𝑌) against house age (𝑋1). From the plot above, we already suspect this may not be ideal, but let us use it as a starting point. lm_house_age_1 =ols('price ~ house_age', data=re2).fit() print(lm_house_age_1.summary()) OLS Regression Results ============================================================================== Dep. Variable: price R-squared: 0.044 Model: OLS Adj. R-squared: 0.042 Method: Least Squares F-statistic: 19.11 Date: Fri, 10 Oct 2025 Prob (F-statistic): 1.56e-05 Time: 09:48:01 Log-Likelihood: -1658.3 No. Observations: 414 AIC: 3321. Df Residuals: 412 BIC: 3329. Df Model: 1 Covariance Type: nonrobust ============================================================================== coef std err t P>|t| [0.025 0.975] ------------------------------------------------------------------------------ Intercept 42.4347 1.211 35.042 0.000 40.054 44.815 house_age -0.2515 0.058 -4.372 0.000 -0.365 -0.138 ============================================================================== Omnibus: 48.404 Durbin-Watson: 1.957 Prob(Omnibus): 0.000 Jarque-Bera (JB): 119.054 104 Skew: 0.589 Prob(JB): 1.40e-26 Kurtosis: 5.348 Cond. No. 39.0 ============================================================================== Notes: [1] Standard Errors assume that the covariance matrix of the errors is correctly specified. From the output, we can tell that the estimated model for Price (𝑌) against Housing age (𝑋1) is: 𝑌 =42.43−0.25𝑋1 The estimates are  𝛽0=42.43and  𝛽1=−0.25. The output includes the 95% confidence intervals for 𝛽0and 𝛽1. The 𝑅2is 0.044, which means that means that only 4.4% of the variation in 𝑌 is explained by 𝑋. This is extremely poor, even though the 𝑝-value for the 𝐹-test is very small (0.000016). A simple interpretation of the model is as follows: For every 1 year increase in house age, there is an average associated decrease in price of 0.25×10,000New Taiwan Dollars. ÁWarning Note that this interpretation has to be taken very cautiously, especially when there are other explanatory variables in the model. Example 5.3 (Price vs. house age estimated line).In linear regression, we almost always wish to use the model to understand what the mean of future observations would be. In this case, we may wish to use the model to understand how the Price changes as house age increases. This is because, based on our formulation, 𝐸(𝑌|𝑋)=𝛽0+𝛽1𝑋 After estimating the parameters, we would have:  𝐸(𝑌|𝑋)=  𝛽0+ 𝛽1𝑋 Thus we can vary the values of 𝑋to study how the mean of 𝑌changes. Here is how we can do so for the model that we have just fit. new_df =sm.add_constant(pd.DataFrame({'house_age' : np.linspace(0,45,100)})) predictions_out =lm_house_age_1.get_prediction(new_df) ax =re2.plot(x='house_age', y='price', kind='scatter', alpha=0.5 ) ax.set_title('Price vs. age'); ax.plot(new_df.house_age, predictions_out.conf_int()[:, 0].reshape(-1), color='blue', linestyle='dashed'); ax.plot(new_df.house_age, predictions_out.conf_int()[:, 1].reshape(-1), color='blue', linestyle='dashed'); ax.plot(new_df.house_age, predictions_out.predicted, color='blue'); 105 Example 5.6 (Price vs. num stores and distance to MRT).For this example, let us work with a reduced model in order to understand how things work. Price will remain the dependent variable, but we shall use distance to MRT (quantitative) and number of nearby convenience stores. However, we shall recode the number of stores as low (or high) corresponding to whether or not there were 4 stores or less (resp. more than 5). re2['num_stores_cat']=['low' if x<= 4else 'high' for xin re2.num_stores] lm_cat_1 =ols('price ~ dist_MRT + num_stores_cat', re2).fit() print(lm_cat_1.summary()) OLS Regression Results ============================================================================== Dep. Variable: price R-squared: 0.502 Model: OLS Adj. R-squared: 0.500 Method: Least Squares F-statistic: 207.3 Date: Fri, 10 Oct 2025 Prob (F-statistic): 5.54e-63 Time: 09:48:02 Log-Likelihood: -1523.3 No. Observations: 414 AIC: 3053. Df Residuals: 411 BIC: 3065. Df Model: 2 Covariance Type: nonrobust ========================================================================================= coef std err t P>|t| [0.025 0.975] ----------------------------------------------------------------------------------------- Intercept 47.8055 0.696 68.679 0.000 46.437 49.174 num_stores_cat[T.low] -7.4086 1.171 -6.325 0.000 -9.711 -5.106 dist_MRT -0.0055 0.000 -11.913 0.000 -0.006 -0.005 ============================================================================== Omnibus: 190.015 Durbin-Watson: 2.138 Prob(Omnibus): 0.000 Jarque-Bera (JB): 2327.960 Skew: 1.618 Prob(JB): 0.00 Kurtosis: 14.157 Cond. No. 4.31e+03 ============================================================================== Notes: [1] Standard Errors assume that the covariance matrix of the errors is correctly specified. [2] The condition number is large, 4.31e+03. This might indicate that there are strong multicollinearity or other numerical problems. The categorical variable has been coded (by Python) as follows: 𝑋4={1, if there were 4 nearby stores or less 0, otherwise As a result, we have estimated two models: • Corresponding to a large number of nearby stores: 𝑌 =47.81−0.0055𝑋2 • Corresponding to a small number of nearby stores: 𝑌 =40.40−0.0055𝑋2 112 5.4.1 Including an Interaction Term A more complex model arises from an interaction between two terms. Here, we shall consider an interaction between a continuous variable and a categorical explanatory variable. Suppose that we have three predictors: height (𝑋1), weight (𝑋2) and gender (𝑋3). As spelt out in the previous section, we should use indicator variables to represent 𝑋3in the model. If we were to include an interaction between gender and weight, we would be allowing for a males and females to have separate coefficients for 𝑋2. Here is what the model would appear as: 𝑌 =𝛽0+𝛽1𝑋1+𝛽2𝑋2+𝛽3𝑋3+𝛽4𝑋2𝑋3+𝑒 Remember that 𝑋3will be 1 for males and 0 for females. The simplified equation for males would be: 𝑌 =(𝛽0+𝛽3)+𝛽1𝑋1+(𝛽2+𝛽4)𝑋2+𝑒 For females, it would be: 𝑌 =𝛽0+𝛽1𝑋1+𝛽2𝑋2+𝑒 Both the intercept and coefficient of 𝑋2are different, for each value of 𝑋3. Recall that in the previous section, only the intercept term was different. Example 5.7 (Interaction between num of stores and distance to MRT).In an earlier model, we fitted the number of stores and distance to MRT to the model. Here, we include an interaction term. lm_cat_2 =ols('price ~ dist_MRT * num_stores_cat', re2).fit() print(lm_cat_2.summary()) OLS Regression Results ============================================================================== Dep. Variable: price R-squared: 0.538 Model: OLS Adj. R-squared: 0.534 Method: Least Squares F-statistic: 158.9 Date: Fri, 10 Oct 2025 Prob (F-statistic): 2.54e-68 Time: 09:48:02 Log-Likelihood: -1508.0 No. Observations: 414 AIC: 3024. Df Residuals: 410 BIC: 3040. Df Model: 3 Covariance Type: nonrobust ================================================================================================== coef std err t P>|t| [0.025 0.975] -------------------------------------------------------------------------------------------------- Intercept 54.8430 1.425 38.499 0.000 52.043 57.643 num_stores_cat[T.low] -14.9553 1.759 -8.504 0.000 -18.412 -11.498 dist_MRT -0.0278 0.004 -6.948 0.000 -0.036 -0.020 dist_MRT:num_stores_cat[T.low] 0.0226 0.004 5.602 0.000 0.015 0.031 ============================================================================== Omnibus: 215.525 Durbin-Watson: 2.117 Prob(Omnibus): 0.000 Jarque-Bera (JB): 3236.645 113 Skew: 1.844 Prob(JB): 0.00 Kurtosis: 16.192 Cond. No. 1.10e+04 ============================================================================== Notes: [1] Standard Errors assume that the covariance matrix of the errors is correctly specified. [2] The condition number is large, 1.1e+04. This might indicate that there are strong multicollinearity or other numerical problems. Notice that we now have the largest adjusted 𝑅2out of all of the models we have fit so far. The model that we have fit consists of two models: 1. For the case that the number of stores is 4 or less: 𝑌 =39.88−0.005𝑋2 2. For the case that the number of stores is more than 4: 𝑌 =54.84−0.028𝑋2. ĹNote Can you interpret the result above based on your intuition or understanding of real estate prices? 5.5 Residual Analysis Recall from earlier that residuals are computed as 𝑟𝑖=𝑌𝑖− 𝑌𝑖 Residual analysis is a standard approach for identifying how we can improve a model. In the case of linear regression, we can use the residuals to assess if the distributional assumptions hold. We can also use residuals to identify influential points that are masking the general trend of other points. Finally, residuals can provided some direction on how to improve the model. 5.5.1 Standardised Residuals It can be shown that the variance of the residuals is in fact not constant! Let us define the hat-matrix as H=X(X′X)−1X′ The diagonal values of Hwill be denoted ℎ𝑖𝑖, for 𝑖=1,…,𝑛. It can then be shown that 𝑉𝑎𝑟(𝑟𝑖)=𝜎2(1−ℎ𝑖𝑖), 𝐶𝑜𝑣(𝑟𝑖,𝑟𝑗)=−𝜎2ℎ𝑖𝑗 As such, we use the standardised residuals when checking if the assumption of Normality has been met. 𝑟𝑖,𝑠𝑡𝑑 =𝑟𝑖 𝜎√1−ℎ𝑖𝑖 If the model fits well, standardised residuals should look similar to a 𝑁(0,1)distribution. In addition, large values of the standardised residual indicate potential outlier points. 114 By the way, ℎ𝑖𝑖 is also referred to as the leverage of a point. It is a measure of the potential influence of a point (on the parameters, and future predictions). ℎ𝑖𝑖 is a value between 0 and 1. For a model with 𝑝parameters, the average ℎ𝑖𝑖 should be should be 𝑝/𝑛. We consider points for whom ℎ𝑖𝑖 >2×𝑝/𝑛to be high leverage points. In the literature and in textbooks, you will see mention of residuals, standardised residuals and studentised residuals. While they differ in definitions slightly, they typically yield the same information. Hence we shall stick to standardised residuals for our course. Example 5.8 (Normality check for lm_age_mrt_1).One of the first checks for Normality is to create a histogram. If the residuals adhere to a Normal distribution, we should observe a symmetric bell-shaped distribution. plt.figure(figsize=(4,3)) r_s =pd.Series(lm_age_mrt_1.resid_pearson) r_s.hist(); 420246 0 50 100 150 200 From above, it appears that the distribution is slightly skewed, and there is one noticeable outlier. A second graphical diagnostic plot that we make is a QQ-plot. A QQ-plot plots the standardized sample quantiles against the theoretical quantiles of a N(0; 1) distribution. If they fall on a straight line, then we would say that there is evidence that the data came from a normal distribution. Especially for unimodal datasets, the points in the middle will fall close to the line. The value of a QQ-plot is in judging if the tails of the data are fatter or thinner than the tails of the Normal. fig =sm.qqplot(r_s[r_s <6], line="q") fig.set_size_inches(4,3); 115 21012 Theoretical Quantiles 4 2 0 2 4 Sample Quantiles Overall, the residuals do indicate a lack of Normal behaviour. Non-normality in the residuals should lead us to view the hypothesis tests with caution. The estimated models are still valid, in the sense that they are optimal. Estimation of the models did not require the assumption of Normality. So far, we have only focused on inspecting the 𝑅2to assess the model quality. 5.5.2 Scatterplots To understand the model fit better, a set of scatterplots are typically made. These are plots of standardised residuals (on the 𝑦-axis) against • fitted values • explanatory variables, one at a time. • potential variables. Residuals are meant to contain only the information that our model cannot explain. Hence, if the model is good, the residuals should only contain random noise. There should be no apparent pattern to them. If we find such a pattern in one of the above plots, we would have some clue as to how we could improve the model. We typically inspect the plots for the following patterns: Figure 5.1: Residuals 116 1. A pattern like the one on the extreme left is ideal. Residuals are randomly distributed around zero; there is no pattern or trend in the plot. 2. The second plot is something rarely seen. It would probably appear if we were to plot residuals against a new variable that is not currently in the model. If we observe this plot, we should then include this variable in the model. 3. This plot indicates we should include a quadratic term in the model. 4. The wedge shape (or funnel shape) indicates that we do not have homoscedascity. The solution to this is either a transformation of the response, or weighted least squares. You will cover these in your linear models class. Example 5.9 (Example: Residual Plots for lm_age_mrt_2).Let us extract and create the residual plots for the second model that we had fit, earlier. plt.figure(figsize=(12,4)); r_s =lm_age_mrt_2.resid_pearson ax=plt.subplot(121) ax.scatter(re2.dist_MRT, r_s, alpha=0.5) ax.set_xlabel('dist_mrt') ax.axhline(y=0, color='red', linestyle='--') ax=plt.subplot(122) ax.scatter(re2.house_age, r_s, alpha=0.5) ax.set_xlabel('house age'); ax.axhline(y=0, color='red', linestyle='--'); 0 1000 2000 3000 4000 5000 6000 dist_mrt 4 2 0 2 4 6 8 0 10 20 30 40 house age 4 2 0 2 4 6 8 While the plot for house age looks acceptable (points are evenly scattered about the red dashed line), the one for distance shows some curvature. This is something we can try to fix, using a transformation of the x-variable. 5.5.3 Influential Points The influence of a point on the inference can be judged by how much the inference changes with and without the point. For instance to assess if point 𝑖is influential on coefficient 𝑗: 1. Estimate the model coefficients with all the data points. 117 2. Leave out the observations (𝑌𝑖,𝑋𝑖)one at a time and re-estimate the model coefficients. 3. Compare the 𝛽’s from step 2 with the original estimate from step 1. While the above method assesses influence on parameter estimates, Cook’s distance performs a similar iteration to assess the influence on the fitted values. Cook’s distance values greater than 1 indicate possibly influential points. There are several ways to deal with influential points. First, we can remove the influential point (or sets of points) and asssess how much the model changes. Based on our understanding of the domain, we can then decide to keep or remove those points. A second approach is to create a dummy variable that identifies those points (individually). Fitting the subsequent model allows all points to be used in estimating standard errors, but quantifies an adjustment for those points. A third approach is to use a robust linear model. This is a model that automatically reduces the influence of aberrant points. This is a good topic to know about - do read up on it if you are keen! Example 5.10 (Example: Influential Points for lm_age_mrt_2).The influence of a point on the inference can be judged by how much the inference changes with and without the point. For instance to assess if point 𝑖is influential on coefficient 𝑗: 1. Estimate the model coefficients with all the data points. 2. Leave out the observations (𝑌𝑖,𝑋𝑖)one at a time and re-estimate the model coefficients. 3. Compare the 𝛽’s from step 2 with the original estimate from step 1. While the above method assesses influence on parameter estimates, Cook’s distance performs a similar iteration to assess the influence on the fitted values. Cook’s distance values greater than 1 indicate possibly influential points. infl =lm_age_mrt_2.get_influence() infl_df =infl.summary_frame() print(infl_df.head()) dfb_Intercept dfb_house_age dfb_x3 dfb_dist_MRT cooks_d \ 0 0.004477 -0.015423 0.004746 0.015011 0.000249 1 -0.004122 0.022294 -0.024375 -0.017434 0.000236 2 0.016558 0.009214 -0.017881 -0.018001 0.000405 3 0.037359 0.020789 -0.040346 -0.040616 0.002054 4 -0.037629 0.024643 -0.013610 0.006609 0.000375 standard_resid hat_diag dffits_internal student_resid dffits 0 -0.350312 0.008050 -0.031558 -0.349937 -0.031524 1 0.319229 0.009195 0.030754 0.318879 0.030720 2 0.638887 0.003955 0.040257 0.638426 0.040228 3 1.438604 0.003955 0.090648 1.440488 0.090767 4 -0.463313 0.006946 -0.038748 -0.462869 -0.038711 118 5.6 Transformation Example 5.11 (Log-transformation).As we observed in the residual plots, the distance-toMRT variable displays a slight curvature. We can fix this by taking a log-transformation of the variable before fitting the model. Below, we include the code to perform this fitting. re2['ldist']=np.log(re2.dist_MRT) lm_age_mrt_3 =ols('price ~ house_age + x3 + num_stores + ldist', data=re2).fit() print(lm_age_mrt_3.summary()) OLS Regression Results ============================================================================== Dep. Variable: price R-squared: 0.597 Model: OLS Adj. R-squared: 0.593 Method: Least Squares F-statistic: 151.6 Date: Fri, 10 Oct 2025 Prob (F-statistic): 2.07e-79 Time: 09:48:02 Log-Likelihood: -1479.4 No. Observations: 414 AIC: 2969. Df Residuals: 409 BIC: 2989. Df Model: 4 Covariance Type: nonrobust ============================================================================== coef std err t P>|t| [0.025 0.975] ------------------------------------------------------------------------------ Intercept 84.0709 4.009 20.971 0.000 76.190 91.951 house_age -0.4941 0.075 -6.578 0.000 -0.642 -0.346 x3 0.8599 0.197 4.355 0.000 0.472 1.248 num_stores 0.7815 0.201 3.886 0.000 0.386 1.177 ldist -6.6595 0.559 -11.910 0.000 -7.759 -5.560 ============================================================================== Omnibus: 233.692 Durbin-Watson: 2.059 Prob(Omnibus): 0.000 Jarque-Bera (JB): 3847.453 Skew: 2.030 Prob(JB): 0.00 Kurtosis: 17.372 Cond. No. 213. ============================================================================== Notes: [1] Standard Errors assume that the covariance matrix of the errors is correctly specified. Now, take some to investigate the following issues: 1. Interpret the coefficient for dist_MRT. 2. Have the issues with the residuals been fixed? 3. What is the difference between how this model uses num_stores, and how lm_cat_1 uses it? 4. How did we choose 25 as the breakpoint for house-age? Is it the ideal one? 119 5.7 Summary, Further topics Linear regression is a very flexible model. It is quick to fit, easily generalisable and much more interpretable than other models. These are some of the reasons why it is still one of the most widely used models in industry. In our short session, we have touched on several practical tips for using regression models. However, take note that regression models can be generalised in many other ways. Here are some models you may want to read up on: • Assuming correlated errors instead of independent error terms • Using splines to include non-linear functions of explanatory variables. • Kernel regression is an even more modern method for including higher-order terms, but at this point we start to lose interpretability • Constrained regression, when we know certain coefficients should be positive, for instance. • Robust linear models to automagically take care of wild outliers. Good reference textbooks for this topic are Draper (1998) and Hastie, Tibshirani, and Friedman (2009). 5.8 References 5.8.1 Website References 1. Taiwan dataset from UCI machine learning repository 2. Stats models documentation 3. Diagnostics 4. On residual plots 120 6 Time Series Analysis 6.1 Exploring Time Series Data One of the first things that we do when we are given a time series is visualise it. The visualisation is meant to give us an indication of what kinds of techniques would be suitable for forecasting it. In this section, we shall learn several methods to visualise a time series dataset. When we visualise a time series, we look out for the following features: 1. Trend: A trend exists when there is a long-term increase or decrease in the data. 2. Level: The level of a series refers to its height on the ordinate axis. 3. Seasonal: A seasonal pattern exists when a series is influenced by factors such as quarters of the year, the month, the day of the week, or time of day. Seasonality is always of a fixed and known period. 4. Cyclic: A cyclic pattern exists when there are rises and falls that are not of a fixed period. import pandas as pd import numpy as np import datetime, calendar from statsmodels.tsa.seasonal import seasonal_decompose,STL from statsmodels.tsa.statespace.tools import diff from statsmodels.tsa.stattools import acf from statsmodels.tsa.forecasting.theta import ThetaModel from statsmodels.tsa.arima.model import ARIMA from statsmodels.tsa.statespace import exponential_smoothing from statsmodels.tsa.api import ( ExponentialSmoothing, SimpleExpSmoothing, Holt, STLForecast ) import pmdarima as pm from scipy.cluster import hierarchy from scipy.spatial.distance import pdist,squareform from ind5003 import ts import matplotlib.pyplot as plt import seaborn as sns Example 6.1 (Basic plots of housing Data).This dataset comes from an R package fma. It contains counts of monthly sales of new one-family houses sold in the USA since 1973. 121 25 50 75 hsales 40 60 Trend 10 0 10 Seasonal 1976 1980 1984 1988 1992 10 0 10 Resid Figure 6.3: Additive decomposition When we inspect a decomposition, there are certain elements that we watch out for. One of them is to compare the range of the 𝑦-axis for the three components. In Figure 6.3, we observe that the rough spread of the trend component (second row) is about 40. On the other hand, the spread of the seasonal component (third row) is roughly 20. This suggests that the seasonal component is a significant driver of the overall variation of the process. The fact that the seasonal pattern does not change over time is an artifact of the decomposition; it is not always a true reflection of the true pattern. Finally, as always, we hope to see the residuals to be centred around 0, with little or no trend, and a constant variability over time. Example 6.6 (Multiplicative decomposition, Aus electricity data).Here is the multiplicative decomposition for the Australian electric data. qau_mult =seasonal_decompose(qau.loc[:, 'kWh'], model='multiplicative', extrapolate_trend='freq') ax =qau_mult.plot(); ax.set_figheight(6) 128 25 50 kWh 25 50 Trend 0.95 1.00 1.05 Seasonal 1960 1970 1980 1990 2000 2010 0.0 0.5 1.0 Resid Take note that in multiplicative decompositions, the residuals are centred around 1, not zero. As you may have noticed, there are some issues with the classical seasonal decomposition algorithms. These include: • an inability to estimate the trend at the ends of the series. • an inability to account for changing seasonal components. • it is not robust to outliers. Example 6.7 (STL decomposition, Aus electricity).Newer algorithms such as STL decomposition, X11 and others attempt to address these issues. They use locally weighted regression models to obtain the trend. These algorithms iterate over the time series several times, so as to ensure that outliers are not affecting the outcome. For full details on STL (Seasonal Trend decomposition with Loess), please refer to the original paper Cleveland et al. (1990). qau_stl =STL(qau.kWh).fit() ax =qau_stl.plot(); ax.set_figheight(6) 129 25 50 kWh 25 50 Trend 0 2 Season 1960 1970 1980 1990 2000 2010 2 0 2 Resid Figure 6.4: STL decomposition, electricity data Notice in Figure 6.4 how the seasonal pattern varies over time. One disadvantage of the STL decomposition, however, is that it can only return an additive decomposition. 6.3 Forecasting 6.3.1 Benchmark methods As in all forecasting methods, it is useful to obtain a baseline forecast before proceeding to more sophisticated techniques. Baseline forecasts are usually obtained from simple, intuitive methods. ĹNote Take note in the sections that follow, the notation 𝑦𝑇+ℎ|𝑇 refers to the forecasted value of the time series at time 𝑇+ℎ, given observations until and including time 𝑇. Here are some baseline methods: 130 A. The simple mean forecast: 𝑦𝑇+ℎ|𝑇 =𝑦1+𝑦2+𝑦3+⋯+𝑦𝑇 𝑇(6.4) B. The naive forecast: 𝑦𝑇+ℎ|𝑇 =𝑦𝑇(6.5) C. The seasonal naive forecast. Suppose that the season contains 𝑑periods (for instance, in monthly data, 𝑑=12). Then the seasonal naive forecast utilises the most recent observed period for the forecast: 𝑦𝑇+ℎ|𝑇 =𝑦𝑇−𝑑+ℎ (6.6) Example 6.8 (Benchmark forecasts housing sales).Suppose we apply some of the above forecasts to the housing sales dataset. We withhold the most recent 2 years of data and forecast those. Here is a plot depicting the forecasts, and the true values. # Set aside the last two years as the test set. #hsales = hsales.drop(columns=['year', 'month']) train_set =hsales.iloc[:-24,] test_set =hsales.iloc[-24:, ] # Obtain the forecast from the training set mean_forecast =ts.meanf(train_set.hsales, 24) snaive_forecast =ts.snaive(train_set.hsales, 24,12) # Plot the predictions and true values ax =train_set.hsales.plot(title='Benchmarks', legend=False, figsize=(12,4.5)) test_set.hsales.plot(ax=ax, legend=False, style='--') mean_forecast.plot(ax=ax, legend=True, style='-') snaive_forecast.plot(ax=ax, legend=False, style='-') plt.legend(labels=['train','test','mean','snaive'], loc='lower right'); 1974 1979 1984 1989 1994 date 30 40 50 60 70 80 90 Benchmarks train test mean snaive Figure 6.5: Benchmark forecasts 131 When assessing forecasts, there are a few different metrics that are typically applied. Each of them has it’s own set of pros and cons. If we denote the predicted value with 𝑦𝑡, then these are the formulas for three of the most common error metrics used in time series forecasting A. RMSE: √ √ √ ⎷1 ℎℎ ∑ 𝑖=1(𝑦𝑡+𝑖− 𝑦𝑡+𝑖)2(6.7) B. MAE 1 ℎℎ ∑ 𝑖=1|𝑦𝑡+𝑖− 𝑦𝑡+𝑖|(6.8) C. Mean Absolute Scaled Error 1 ℎℎ ∑ 𝑖=1 |𝑦𝑡+𝑖− 𝑦𝑡+𝑖| 1 𝑇−1∑𝑇 𝑡=2|𝑦𝑡−𝑦𝑡−1|(6.9) The RMSE and MAE are scale dependent errors. It is difficult to compare the errors across series, or to aggregate errors across different time series with it. Due to the square in the formula, the RMSE is quite sensitive to outliers. The MAE is more robust to outliers. The MASE is a scaled error - it allows us to compare the forecasting performance across time series. The other two metrics depend on the scale of the time series and hence do not allow us to make such comparisons. Example 6.9 (Benchmark forecast errors).Let us assess the forecasts made in the previous example, using RMSE and MAE. for xin [ts.rmse, ts.mae]: print(f"{x.__name__},mean: {x(test_set.hsales.values, mean_forecast.values):.3f}") print(f"{x.__name__},snaive: {x(test_set.hsales.values, snaive_forecast.values):.3f}") print('---') rmse,mean: 9.023 rmse,snaive: 5.906 --- mae,mean: 7.562 mae,snaive: 4.792 --- In our simple example, the seasonal naive model outperforms the simple mean forecast according to all the metrics but in general, things are not always this clear-cut. For reference, the MASE for the seasonal naive forecast is as follows. ts.mase(test_set.hsales.values, snaive_forecast.values, train_set.values, seasonality=12) 1.5154940449933834 132 6.3.2 ARIMA Models Now let us turn to a huge class of models that have been utilised in time series forecasting since the 1960’s. They are known as ARIMA models. ARIMA stands for AutoRegressive Integrated Moving Average models. These models are appropriate for stationary processes. Stationarity is a technical term that refers to processes • that have a constant mean. This means that the process merely fluctuates about a fixed level over time. • whose covariance function does not change over time. This means that, for a fixed ℎ, the covariance between 𝑦𝑡and 𝑦𝑡+ℎ is the same for all 𝑡. • whose variance is constant over time. How can we tell if a process is stationary or not? The following time series is not stationary. Why? Example 6.10 (Dow Jones index).Consider the following time series of the Dow Jones index. dj =pd.read_csv('data/dj.csv') dj.plot(legend=False, title='Dow Jones Index', figsize=(12,4)); 0 50 100 150 200 250 300 3600 3700 3800 3900 4000 Dow Jones Index It is nonstationary because the mean (level) of the time series is not constant. However, the following differenced version of the same series is: Δ𝑦𝑡=𝑦𝑡−𝑦𝑡−1 (6.10) ddj =diff(dj) ddj.plot(legend=False, title='Differenced Dow Jones', figsize=(12,4)); 133 0 50 100 150 200 250 300 100 75 50 25 0 25 50 75 Differenced Dow Jones ARIMA models revolve around the idea that if we have a non-stationary series, we can transform it into a stationary one with a suitable number of differencing. A more general way of studying if a series is stationary is to plot it’s AutoCorrelation Function (ACF). The ARIMA method directly models the ACF. That is why it is so important for this class of models. The ACF of a stationary process should “die down” quickly. Figure 6.6 shows the ACF of the Dow Jones data, before and after differencing. plt.figure(figsize=(8,6)) plt.subplot(211) plt.stem(acf(dj, fft=False)) plt.title("Non-differenced") plt.subplot(212) plt.stem(acf(ddj, fft=False)); plt.title("Differenced Series"); 134 0 5 10 15 20 25 0.0 0.2 0.4 0.6 0.8 1.0 Non-differenced 0 5 10 15 20 25 0.0 0.2 0.4 0.6 0.8 1.0 Differenced Series Figure 6.6: ACF, before and after differencing Now suppose that, starting from our original series 𝑦𝑡, we difference it a sufficient number of times and obtain a stationary series. Let’s call this 𝑦′𝑡. The ARIMA model assumes that 𝑦′𝑡=𝑐+𝜙1𝑦′𝑡−1+𝜙2𝑦′𝑡−2+⋯+𝜙𝑝𝑦′𝑡−𝑝+𝜃1𝑒𝑡−1+⋯+𝜃𝑞𝑒𝑡−𝑞 (6.11) • The 𝑒𝑗correspond to unobserved innovations. They are typically assumed to be independent across time with a common variance. • The 𝜙and 𝜃terms are unknown coefficients to be estimated. • If 𝑦′𝑡was obtained by performing 𝑑successive differencings, then the above ARIMA model is referred to as an ARIMA(𝑝,𝑑,𝑞)model. In olden days, the 𝑝,𝑑and 𝑞parameters were picked by the analyst after inspecting the ACF, PACF, and time plots of the differenced series. Today, we can iterate through a large number of them and pick the best according to a well-established criteria (AIC). Example 6.11 (Auto ARIMA on Aus electricity).Let us try out the ARIMA auto-fitting algorithm on the qau electricity usage dataset. We shall set aside the last three years of data as the test set. Recall that this is quarterly data. First, we establish the seasonal naive benchmark error for this dataset. 135 train2 =qau.kWh[:-12] test2 =qau.kWh[-12:] snaive_f =ts.snaive(train2, 12,4) print(f"The RMSE is approximately {ts.rmse(test2.values, snaive_f.values):.3f}.") The RMSE is approximately 2.545. Let us see if the automatically fitted ARIMA models can do better than this. arima_m1 =pm.auto_arima(train2.values, seasonal=True, m=4, test='adf', trace=False, suppress_warnings=True) Here is the summary of the final chosen model. print(arima_m1.summary()) SARIMAX Results ============================================================================================== Dep. Variable: y No. Observations: 206 Model: SARIMAX(1, 1, 1)x(2, 0, [1, 2], 4) Log Likelihood -174.692 Date: Tue, 14 Oct 2025 AIC 365.384 Time: 14:10:30 BIC 391.968 Sample: 0 HQIC 376.136 - 206 Covariance Type: opg ============================================================================== coef std err z P>|z| [0.025 0.975] ------------------------------------------------------------------------------ intercept 0.0023 0.004 0.625 0.532 -0.005 0.009 ar.L1 0.5358 0.081 6.621 0.000 0.377 0.694 ma.L1 -0.8998 0.053 -17.132 0.000 -1.003 -0.797 ar.S.L4 0.0945 0.215 0.439 0.661 -0.327 0.516 ar.S.L8 0.8794 0.212 4.146 0.000 0.464 1.295 ma.S.L4 0.3137 0.208 1.510 0.131 -0.093 0.721 ma.S.L8 -0.5930 0.114 -5.220 0.000 -0.816 -0.370 sigma2 0.3182 0.028 11.178 0.000 0.262 0.374 =================================================================================== Ljung-Box (L1) (Q): 0.03 Jarque-Bera (JB): 25.07 Prob(Q): 0.87 Prob(JB): 0.00 Heteroskedasticity (H): 9.13 Skew: -0.13 Prob(H) (two-sided): 0.00 Kurtosis: 4.69 =================================================================================== Warnings: [1] Covariance matrix calculated using the outer product of gradients (complex-step). 136 The output indicates that a seasonal ARIMA model has been found to be the best one. The column coef contains estimates of the parameters in the model. The section below contains information on statistical tests • for Normality on the residuals (Ljung-Box and Jarque-Bera), and • for constant variance. In addition to tests, we should always create plots to assess the validity of residuals. Once we have found the best-fitting model, we do what we do with any other model in data analytics: we have to inspect the residuals. The residuals should look like trash to us - there should be no clues about the data in them. With reference to Figure 6.7, the standardized residual in the top left should be consistently wide; it isn’t. In fact this information is also reflected in the model output summary (see Heteroskedasticity (H)). This suggests some sort of transformation of the data before modeling might be appropriate. The two plots on the off-diagonal are meant for us to assess if the residuals are Normally distributed with mean 0. They do indeed look like it. The final plot, in the bottom right, displays an ACF with no spikes, indicating that the residuals are uncorrelated. This is precisely what we wished to see. arima_m1.plot_diagnostics(figsize=(12,6)); 0 25 50 75 100 125 150 175 200 2 0 2 4Standardized residual 3210123 0.0 0.1 0.2 0.3 0.4 Histogram plus estimated density Hist KDE N(0,1) 21012 Theoretical Quantiles 2 0 2 4 Sample Quantiles Normal Q-Q 0246810 1.0 0.5 0.0 0.5 1.0 Correlogram Figure 6.7: ARIMA residual diagnostics Finally, we assess the error on the test set. The out-of-sample performance is better than the naive methods. 137 Method: OLS/SES Deseasonalized: True Date: Tue, 14 Oct 2025 Deseas. Method: Multiplicative Time: 14:10:33 Period: 12 Sample: 01-01-1978 - 12-01-2019 Parameter Estimates ======================== Parameters ------------------------ b0 2628.8595906871597 alpha 0.8383848338839772 ------------------------ Here is what the two components from the model look like, in comparison to the seasonally adjusted series (based on a multiplicative decomposition). ĹNote The calculations in the next cell are based on formulas for the theta decomposition that we won’t get into. They are just for illustration of how the decomposition is in terms of long and short term trends, rather than the cross-sectional decompositions we have been dealing with so far. The orange dashed line in Figure 6.12 denotes the long-term trend that the model has identified, while the green dashed line denotes the short-term trend. The blue solid line is the seasonally adjusted time series. # carry out mutliplicative decomposition in order to obtain seasonally adjusted series sea2_mult =seasonal_decompose(sea2, model='multiplicative', extrapolate_trend='freq') # get long-term trend slope and intercept a0 =(sea2.mean() -theta_mod_fit.params['b0']*(503)/2)/1e6 xvals =np.arange(1,505) yvals1 =a0 +theta_mod_fit.params['b0']*(xvals -1) yvals1 =pd.Series(data=yvals1, index=sea2.index) # get short term trend series a2 = -1*a0 b2 = -1*theta_mod_fit.params['b0'] yvals2 =a2 +b2*(xvals -1)+2*(sea2_mult.resid +sea2_mult.trend) yvals2 =pd.Series(data=yvals2, index=sea2.index) fig =(sea2_mult.resid +sea2_mult.trend).plot(label='Seasonally adjusted', legend=True, figsize=(12,4)); yvals1.plot(legend=True, label='Long-term trend', style="--") yvals2.plot(legend=True, label='Short-term trend', style="--"); fig.set_title("Theta decomposition"); 144 1979 1984 1989 1994 1999 2004 2009 2014 2019 Date 0.00 0.25 0.50 0.75 1.00 1.25 1.50 1.75 1e6 Theta decomposition Seasonally adjusted Long-term trend Short-term trend Figure 6.12: Theta decomposition Finally, we visualise the forecasts from this theta decomposition model. We only plot the final 360 observations for the visualisation. theta_forecasts =pd.DataFrame( { "original": sea2, "theta forecast": theta_mod_fit.forecast(24) } ) theta_forecasts.tail(360).plot(figsize=(12,4)); 1994 1999 2004 2009 2014 2019 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 1e6 original theta forecast 6.3.5 Forecasting with Seasonal Decomposition Continuing with the tourism data, we demonstrate how we can utilise a seasonal decomposition (not theta decomposition) model to make forecasts. Recall that a seasonal decomposition breaks the series down into the trend, seasonal and residual components. Seasonal decompositions are typically used to study a time series, but they can also be used to make forecasts. One approach is to use either ARIMA or ETS to predict the seasonally adjusted component, and then to use 145 a seasonal naive model to predict the seasonal component. These two forecasts will then be combined to create the final forecast. Example 6.14 (Tourism forecasts with decomposition).Here is the code that will use and STL model to decompose the series, and fit a pre-chosen ARIMA(3,1,2) model to the seasonally adjusted series. stlf =STLForecast(sea2, ARIMA, model_kwargs={"order": (3,1,2)}) res =stlf.fit( ) forecasts2 =pd.DataFrame( { "original": sea2, "dcmp forecast": res.forecast(24) } ) forecasts2.tail(360).plot(figsize=(12,4)); C:\Users\stavg\Documents\courses\ind5003-book\env\lib\site-packages\statsmodels\tsa\statespace\sarimax.py:966: UserWarning: Non-stationary starting autoregressive parameters found. Using zeros as starting parameters. 1994 1999 2004 2009 2014 2019 0.25 0.50 0.75 1.00 1.25 1.50 1.75 1e6 original dcmp forecast 6.4 Miscellaneous Topics 6.4.1 Time Series Clustering Example 6.15 (Time series clustering, US employment data).The Forecasting: Principles and Practice contains time series data on employment in various sectors in the US. There are 148 unique series in the dataset. us_employment =pd.read_csv("data/us_employment.csv", parse_dates=[0], date_format="%Y %b") series_ids =us_employment.Series_ID.unique() us_employment.sample(n=8) 146 Month Series_ID Title Employed 66865 1939-05-01 CEU4348300001 Transportation and Warehousing: Water Transpor... NaN 100457 1993-03-01 CEU6054160001 Professional and Business Services: Management... 376.0 44668 1946-11-01 CEU4000000001 Trade, Transportation, and Utilities 9411.0 70249 1979-02-01 CEU4348600001 Transportation and Warehousing: Pipeline Trans... NaN 141495 1940-10-01 PAYNSA All Employees, Total Nonfarm 33806.0 8896 1953-08-01 CEU1021300001 Mining and Logging: Support Activities for Mining NaN 7779 1941-04-01 CEU1021210001 Mining and Logging: Coal Mining NaN 120186 1941-07-01 CEU7000000001 Leisure and Hospitality 2121.0 Here is a plot of one of the time series in the dataset. us_employment[us_employment.Title == "Total Private"].plot(x='Month', y='Employed', figsize=(12,4)); 1939 1949 1959 1969 1979 1989 1999 2009 2019 Month 40000 60000 80000 100000 120000 Employed Since we are going to perform clustering on the time series data, our first step is to remove the missing values. us_full =us_employment[(us_employment.Month >= datetime.datetime(2002,9,30)) &(us_employment.Month <datetime.datetime(2018,1,1) ) ] us2 =us_full.pivot(index='Series_ID', columns='Month', values="Employed") us2_array =us2.to_numpy() Now we can begin the clustering procedure. The first line in the code below generates a 148 x 148 distance matrix. out =pdist(us2_array, metric='correlation') lm1 =hierarchy.linkage(out, method='average', optimal_ordering=True) plt.figure(figsize=(12,5)) hierarchy.dendrogram(lm1, p=3, truncate_mode='level',color_threshold=True); 147 68 7 8 (41) (18) 42 77 137 (74) (3) 100 (5) 0.0 0.2 0.4 0.6 0.8 1.0 1.2 We reorder the columns and rows in the matrix so that similar time series appear next to one another. X_ord =us2_array[hierarchy.leaves_list(lm1)] corr_mat_ord =np.corrcoef(X_ord) plt.figure(figsize=(15,15)) sns.heatmap(corr_mat_ord, vmin=-1, vmax=1, cmap='coolwarm_r', center=0); 148 0 3 6 9 12 15 18 21 24 27 30 33 36 39 42 45 48 51 54 57 60 63 66 69 72 75 78 81 84 87 90 93 96 99 102 105 108 111 114 117 120 123 126 129 132 135 138 141 144 147 0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40 42 44 46 48 50 52 54 56 58 60 62 64 66 68 70 72 74 76 78 80 82 84 86 88 90 92 94 96 98 100 102 104 106 108 110 112 114 116 118 120 122 124 126 128 130 132 134 136 138 140 142 144 146 1.00 0.75 0.50 0.25 0.00 0.25 0.50 0.75 1.00 Figure 6.13: Heatmap of correlation matrix In the next few plots, we use the heatmap in Figure 6.13 to pick out series that are close to one another. Series in blue segments of the matrix are highly correlated; here we inspect them to try to characterise the similarity. The series in Figure 6.14 seem to have a dip around the period 2009 - 2011, while those in Figure 6.15 have a cyclical pattern. us2_series =us2.index.to_list() ss =[us2_series[x] for xin [22,24,32,33,34]] us2.T.loc[:, ss].plot(figsize=(12,4)); plt.legend(loc='upper right'); 149 2003 2005 2007 2009 2011 2013 2015 2017 Month 400 600 800 1000 1200 CEU3133100001 CEU3133300001 CEU3133600101 CEU3133700001 CEU3133900001 Figure 6.14: Series from top left corner of heatmap ss =[us2_series[x] for xin [143,144,145]] us2.T.loc[:, ss].plot(figsize=(12,4)); 2003 2005 2007 2009 2011 2013 2015 2017 Month 6000 8000 10000 12000 14000 Series_ID CEU9093000001 CEU9093161101 CEU9093200001 Figure 6.15: Series from bottom right corner of heatmap 6.5 Summary We have briefly covered the following time series concepts: Exploratory data analysis, assessing point forecasts and a few commonly used models. To finish up the section, we also applied clustering to time series. With today’s code, the model-fitting routines are easier to run. As always, the onus is on us as analysts to inspect the residuals from the models to understand the time series we are working with. 6.6 References 6.6.1 Stats models pages 1. Main page 2. Forecasting with statsmodels 150 3. Theta model 6.6.2 Forecasting principles and practice Although the code for this textboook uses R, the concepts are very well explained. It is written by one of the foremost experts in time series forecasting methods. The reference is Hyndman and Athanasopoulos (2018). •Forecasting: Principles and Practice 151 7 Simulation 7.1 Random Variables In statistics, we use random variables to describe the probabilistic behaviour of phenomenon. Random variables are real numbers that represent the phenomena we observe. For instance, we might let 𝑋represent a coin toss, with 𝑋=1representing Heads and 𝑋=0representing Tails. Random variables come with a rule (or function) that prescribes the probabilities with which it takes on particular values. For instance, if we had a fair coin, then the rule would be that 𝑃(𝑋=1)=𝑃(𝑋=0)=1 2 On the other hand, a biased coin might follow the rule 𝑃(𝑋=1)=2 3, 𝑃(𝑋=0)=1 3 This rule tells us which events are more likely, and which are less likely. from math import pi import pprint import mesa import seaborn as sns import matplotlib.pyplot as plt import nltk from nltk import word_tokenize from numpy.random import default_rng import numpy as np import pandas as pd from scipy import stats from scipy.stats import binom, bernoulli, norm, expon, uniform 7.1.1 Discrete Random Variables Discrete random variables take on only a countable number of values. Examples are: • A coin toss ({0,1}) 152 • The number of taxis passing by a particular junction between 12noon and 1pm. ({0,1,2,…,}) • The number of coin tosses until we observe Heads. ({1,2,3,…,}) Discrete random variables are defined by their probability mass function (pmf), which is just a table or a function describing 𝑃(𝑋 =𝑖)for all possible 𝑖values. Once we know the pmf of a random variable, we know everything about it - the mean, variance, quantiles, maximum values, etc. We can visualise a pmf using a bar-chart. Here is the pmf for the above biased coin. plt.figure(figsize=(5,3)) plt.bar([0,1], [1/3,2/3], tick_label=['0','1'], width=0.3); plt.xlim(-1,2); plt.ylim(0,1); plt.title('pmf of \'X\' random variable' ); 0 1 0.0 0.2 0.4 0.6 0.8 1.0 pmf of 'X' random variable Here is the pmf for a random variable representing the total number of Heads after 10 tosses of that same coin. Suppose we call that new random variable 𝑌. plt.figure(figsize=(5,3)) probs =binom.pmf(np.arange(0,11), n=10, p=2/3) plt.bar(np.arange(0,11), probs);plt.ylim(0,1); plt.title('pmf of no. of coin tosses'); 153