Line of Best Fit - Portfolio

Line of Best Fit #

Creating a line of best fit to predict future country population, optimizing line hyperparameters through a genetic algorithm and error function.

.rmd knitted with execution results available for download here

Part 1: Download the file#

1
# url = "https://population.un.org/wpp/Download/Files/1_Indicators%20(Standard)/CSV_FILES/WPP2022_Demographic_Indicators_Medium.zip"
2

3
# download.file(url, "output.zip")

Part 2: Unzip the file, put into dataframe#

1
# unzip
2
# zip = "output.zip"
3
# unzip(zip)
4

5
# put into dataframe
6
# file_name = "WPP2022_Demographic_Indicators_Medium.csv"
7
# full_data = read.csv(file_name)
8

9
full_data = read.csv(url("https://raw.githubusercontent.com/julien-arino/math-of-data-science/main/CODE/WPP2022_Demographic_Indicators_Medium.csv"))
10

11
# filter out irrelevant data
12
data = full_data[full_data$Location == "Canada" & full_data$Time <= 2023,]

Part 3: Plot the total population of that country#

Multiply the population by 1000 to get the real number, not in a measure of thousands

1
# sort data by time ascending
2
sorted = data[order(data$Time),]
3

4
# plot data
5
plot(x=data$Time ,y=data$TPopulation1Jan*1000, xlab="Year", ylab="Population")

Part 4a: Fit lines to the plot#

Start with fitting the line: y = a + bx#

Initially, create the matrix ‘A’, which consists of two columns of the coefficients for ‘a’ and ‘b’

The coefficients for the ‘a’ value is just 1’s, so fill the first column with a 1 for every ‘x’ value in the plot

The coefficients for the ‘b’ value are the ‘x’ values, so the second column is all of the values in ‘data$Time’

1
A = matrix(c(rep(1, length(data$Time)), data$Time), nr=length(data$Time), nc=2)
2
A

1
##       [,1] [,2]
2
##  [1,]    1 1950
3
##  [2,]    1 1951
4
##  [3,]    1 1952
5
##  [4,]    1 1953
6
##  [5,]    1 1954
7
##  [6,]    1 1955
8
##  [7,]    1 1956
9
##  [8,]    1 1957
10
##  [9,]    1 1958
11
## [10,]    1 1959
12
## [11,]    1 1960
13
## [12,]    1 1961
14
## [13,]    1 1962
15
## [14,]    1 1963
16
## [15,]    1 1964
17
## [16,]    1 1965
18
## [17,]    1 1966
19
## [18,]    1 1967
20
## [19,]    1 1968
21
## [20,]    1 1969
22
## [21,]    1 1970
23
## [22,]    1 1971
24
## [23,]    1 1972
25
## [24,]    1 1973
26
## [25,]    1 1974
27
## [26,]    1 1975
28
## [27,]    1 1976
29
## [28,]    1 1977
30
## [29,]    1 1978
31
## [30,]    1 1979
32
## [31,]    1 1980
33
## [32,]    1 1981
34
## [33,]    1 1982
35
## [34,]    1 1983
36
## [35,]    1 1984
37
## [36,]    1 1985
38
## [37,]    1 1986
39
## [38,]    1 1987
40
## [39,]    1 1988
41
## [40,]    1 1989
42
## [41,]    1 1990
43
## [42,]    1 1991
44
## [43,]    1 1992
45
## [44,]    1 1993
46
## [45,]    1 1994
47
## [46,]    1 1995
48
## [47,]    1 1996
49
## [48,]    1 1997
50
## [49,]    1 1998
51
## [50,]    1 1999
52
## [51,]    1 2000
53
## [52,]    1 2001
54
## [53,]    1 2002
55
## [54,]    1 2003
56
## [55,]    1 2004
57
## [56,]    1 2005
58
## [57,]    1 2006
59
## [58,]    1 2007
60
## [59,]    1 2008
61
## [60,]    1 2009
62
## [61,]    1 2010
63
## [62,]    1 2011
64
## [63,]    1 2012
65
## [64,]    1 2013
66
## [65,]    1 2014
67
## [66,]    1 2015
68
## [67,]    1 2016
69
## [68,]    1 2017
70
## [69,]    1 2018
71
## [70,]    1 2019
72
## [71,]    1 2020
73
## [72,]    1 2021
74
## [73,]    1 2022
75
## [74,]    1 2023

Now make the vector ‘b’, which is just the ‘y’ values of the plot

1
b = data$TPopulation1Jan

Now create the matrix multiplication of ‘A’ transpose times ‘A’, and get it’s inverse

1
AtA = t(A) %*% A
2

3
# Look at the result, and it's determinate
4
AtA

1
##        [,1]      [,2]
2
## [1,]     74    147001
3
## [2,] 147001 292051249

1
det(AtA)

1
## [1] 2498425

1
# Invert the result
2
invAtA = solve(AtA)

Since the determinate of ‘AtA’ is non-zero, and the columns that make up ‘A’ are clearly linearly independent, there is a unique solution.

The unique solution is the inverse of ‘AtA’, times the transpose of ‘A’, times ‘b’.

1
solution_linear = invAtA %*% t(A) %*% b
2
solution_linear

1
##              [,1]
2
## [1,] -627577.8531
3
## [2,]     329.1654

Now fit the line: y = a + bx + cx^2#

It’s pretty much the same, only ‘A’ has a third column consisting of the coefficients for ‘c’, namely a column of ‘x^2’

However, due to the large values of ‘x^2’, that will continue to get larger as we work towards a solution for our quadratic line, we could run into problems when we try and find the inverse, namely the value of our determinate will get so extreme that the result will be functionally 0.

This can be solved by modifying the values of ‘x’, so long as they maintain the gaps between each value there will be no problem.

So by subtracting the minimum ‘x’ value from all ‘x’ values, we get smaller numbers that maintain the gap between the original ‘x’ values, thereby solving the issue of the extreme determinate.

1
A2 = matrix(c(rep(1, length(data$Time)), data$Time-min(data$Time), (data$Time-min(data$Time))^2), nr=length(data$Time), nc=3)
2
A2

1
##       [,1] [,2] [,3]
2
##  [1,]    1    0    0
3
##  [2,]    1    1    1
4
##  [3,]    1    2    4
5
##  [4,]    1    3    9
6
##  [5,]    1    4   16
7
##  [6,]    1    5   25
8
##  [7,]    1    6   36
9
##  [8,]    1    7   49
10
##  [9,]    1    8   64
11
## [10,]    1    9   81
12
## [11,]    1   10  100
13
## [12,]    1   11  121
14
## [13,]    1   12  144
15
## [14,]    1   13  169
16
## [15,]    1   14  196
17
## [16,]    1   15  225
18
## [17,]    1   16  256
19
## [18,]    1   17  289
20
## [19,]    1   18  324
21
## [20,]    1   19  361
22
## [21,]    1   20  400
23
## [22,]    1   21  441
24
## [23,]    1   22  484
25
## [24,]    1   23  529
26
## [25,]    1   24  576
27
## [26,]    1   25  625
28
## [27,]    1   26  676
29
## [28,]    1   27  729
30
## [29,]    1   28  784
31
## [30,]    1   29  841
32
## [31,]    1   30  900
33
## [32,]    1   31  961
34
## [33,]    1   32 1024
35
## [34,]    1   33 1089
36
## [35,]    1   34 1156
37
## [36,]    1   35 1225
38
## [37,]    1   36 1296
39
## [38,]    1   37 1369
40
## [39,]    1   38 1444
41
## [40,]    1   39 1521
42
## [41,]    1   40 1600
43
## [42,]    1   41 1681
44
## [43,]    1   42 1764
45
## [44,]    1   43 1849
46
## [45,]    1   44 1936
47
## [46,]    1   45 2025
48
## [47,]    1   46 2116
49
## [48,]    1   47 2209
50
## [49,]    1   48 2304
51
## [50,]    1   49 2401
52
## [51,]    1   50 2500
53
## [52,]    1   51 2601
54
## [53,]    1   52 2704
55
## [54,]    1   53 2809
56
## [55,]    1   54 2916
57
## [56,]    1   55 3025
58
## [57,]    1   56 3136
59
## [58,]    1   57 3249
60
## [59,]    1   58 3364
61
## [60,]    1   59 3481
62
## [61,]    1   60 3600
63
## [62,]    1   61 3721
64
## [63,]    1   62 3844
65
## [64,]    1   63 3969
66
## [65,]    1   64 4096
67
## [66,]    1   65 4225
68
## [67,]    1   66 4356
69
## [68,]    1   67 4489
70
## [69,]    1   68 4624
71
## [70,]    1   69 4761
72
## [71,]    1   70 4900
73
## [72,]    1   71 5041
74
## [73,]    1   72 5184
75
## [74,]    1   73 5329

The vector ‘b’ remains unchanged

1
AtA2 = t(A2) %*% A2
2

3
# Look at the result, and it's determinate
4
AtA2

1
##        [,1]    [,2]      [,3]
2
## [1,]     74    2701    132349
3
## [2,]   2701  132349   7295401
4
## [3,] 132349 7295401 428943109

1
det(AtA2)

1
## [1] 3.0772e+13

1
# Invert the result
2
invAtA2 = solve(AtA2)
3
invAtA2

1
##               [,1]          [,2]          [,3]
2
## [1,]  1.152774e-01 -6.273115e-03  7.112376e-05
3
## [2,] -6.273115e-03  4.622882e-04 -5.926980e-06
4
## [3,]  7.112376e-05 -5.926980e-06  8.119150e-08

Now compute the solution which is the inverse of ‘AtA2’, times the transpose of ‘A2’, times ‘b’

1
solution_quadratic = invAtA2 %*% t(A2) %*% b
2
solution_quadratic

1
##               [,1]
2
## [1,] 14104.2697665
3
## [2,]   345.0353083
4
## [3,]    -0.2173958

Part 4b: Fit lines to the plot using Least Squares Fit#

Starting with y = a + bx#

Start by creating a function to model a line.

Then we need to find the error between this line and the actual values.

This is done by taking the difference between the actual points and the modeled points

1
# Create a function to output the estimated values
2
# y = a + bx
3
linear = function(x, a, b){
4
  return(a + b*x)
5
}
6

7
# Create the error function for a linear function
8
# e_i = y_i - y_i~
9
# total is sqrt of sum of e_i ^ 2
10
error_linear = function(a, b, points){
11
  yt = linear(points$x, a, b)
12
  e = points$y - yt
13
  return(sqrt(sum(e^2)))
14
}

But we can’t actually do anything with this unless we have the parameters needed to model a line.

To solve this we can use a genetic algorithm to automatically optimize our hyper parameters using our error function.

1
# Create points list for model
2
points = list()
3
points$x = c(data$Time-min(data$Time))
4
points$y = c(data$TPopulation1Jan)*1000
5

6
# installs 'GA' package if not already installed
7
if(!require(GA)){
8
  install.packages("GA")
9
  library(GA)
10
}

1
## Loading required package: GA

1
## Loading required package: foreach

1
## Loading required package: iterators

1
## Package 'GA' version 3.2.3
2
## Type 'citation("GA")' for citing this R package in publications.

1
##
2
## Attaching package: 'GA'

1
## The following object is masked from 'package:utils':
2
##
3
##     de

1
library(GA)
2
GA = ga(type = "real-valued",
3
        fitness = function(val) {return(-error_linear(a = val[1], b = val[2], points))},
4
        suggestions = c(a = 1.4e+07, b = 300000),
5
        lower = c(0, 0), upper = c(4e+07, 1.5e+10),
6
        popSize = 200, maxiter = 500)
7
plot(GA)

And the values for ‘a’ and ‘b’ are:

1
GA@solution

1
##            x1       x2
2
## [1,] 13629086 347462.6

1
-GA@fitnessValue

1
## [1] 4271004

Now for y = a + bx + cx^2#

Simply create a new quadratic function and an error function to go along with it.

The formula for the error stays the same.

1
# Create a function to output the estimated values
2
# y = a + bx + cx^2
3
quadratic = function(x, a, b, c){
4
  return(a + b*x + c*(x^2))
5
}
6

7
# Create the error function for a quadratic function
8
# e_i = y_i - y_i~
9
# total is sqrt of sum of e_i ^ 2
10
error_quadratic = function(a, b, c, points){
11
  yt = quadratic(points$x, a, b, c)
12
  e = points$y - yt
13
  return(sqrt(sum(e^2)))
14
}

Then just preform the genetic algorithm with the new error function

1
GA2 = ga(type = "real-valued",
2
        fitness = function(val) {return(-error_quadratic(a = val[1], b = val[2], c=val[3], points))},
3
        suggestions = c(a = 1.4e+07, b = 3000, c = 10),
4
        lower = c(0, 0, 0), upper = c(2e+07, 1.5e+7, 3000),
5
        popSize = 200, maxiter = 500)
6
plot(GA2)

1
GA2@solution

1
##            x1       x2       x3
2
## [1,] 13736523 339490.7 8.124464

1
-GA2@fitnessValue

1
## [1] 3614614

Part 5: Plot lines on graph#

1
# use the modified x-values
2
x = data$Time-min(data$Time)
3

4
# plot data every time
5

6
# plot y = a + bx
7
plot(x=data$Time ,y=data$TPopulation1Jan*1000, xlab="Year", ylab="Population")
8
title("y = a + bx")
9
lines(x=data$Time, y=1000*(solution_linear[1] + solution_linear[2]*data$Time), type="l", lwd=2, col="red")

1
# plot y = a + bx + cx^2
2
plot(x=data$Time ,y=data$TPopulation1Jan*1000, xlab="Year", ylab="Population")
3
title("y = a + bx + cx^2")
4
lines(x=data$Time, y=1000*(solution_quadratic[1] + solution_quadratic[2]*x + solution_quadratic[3]*(x ^ 2)), type="l", lwd=1, col="blue")

1
# plot y = a + bx least squares
2
plot(x=data$Time ,y=data$TPopulation1Jan*1000, xlab="Year", ylab="Population")
3
title("y = a + bx least squares")
4
lines(x=data$Time, y=(linear((points$x), GA@solution[1], GA@solution[2])), lwd=2, col="black")

1
# plot y = a + bx + cx^2 least squares
2
plot(x=data$Time ,y=data$TPopulation1Jan*1000, xlab="Year", ylab="Population")
3
title("y = a + bx + cx^2 least squares")
4
lines(x=data$Time, y=(quadratic((points$x), GA2@solution[1], GA2@solution[2], GA2@solution[3])), lwd=2, col="green")

Part 6: Predict the population at 2030#

Predicting the population at 2030 just involves plugging in 2030 for the ‘x’ values used in the ‘y=’ part of the line plots.

Though, you have to remember to shift 2030 by the earliest year value for the quadratic prediction and the least squares prediction.

1
linear_pred = 1000*(solution_linear[1] + solution_linear[2]*2030)
2
quad_pred = 1000*(solution_quadratic[1] + solution_quadratic[2]*(2030-min(data$Time)) + solution_quadratic[3]*((2030-min(data$Time)) ^ 2))
3
least_lpred = linear((2030-min(data$Time)), GA@solution[1], GA@solution[2])
4
least_qpred = quadratic((2030-min(data$Time)), GA2@solution[1], GA2@solution[2], GA2@solution[3])
5

6
linear_pred

1
## [1] 40627942

1
quad_pred

1
## [1] 40315761

1
least_lpred

1
## [1] 41426094

1
least_qpred

1
## [1] 40947774

It should be noted that the original data has predictions for Canada’s future population up to the year 2100, so we can compare our predicted result, to the predicted result in the .csv.

1
pred = full_data[full_data$Location == "Canada" & full_data$Time == 2030,]
2
pred$TPopulation1Jan * 1000

1
## [1] 40852177

From this we can see that the closest model to the given prediction is the least squares quadratic model.

Line of Best Fit#