#set document(title: "6.2 Classification Using Machine Learning", author: "OpenStax / XYZ Homework")
#set page(width: 8.5in, height: auto, margin: 1in)
#import "@preview/cetz:0.5.2"
#set text(font: ("STIX Two Text", "Libertinus Serif", "New Computer Modern"), size: 10.5pt, lang: "en")
#show math.equation: set text(font: ("STIX Two Math", "New Computer Modern Math"))
#set par(justify: true, leading: 0.62em, spacing: 0.9em)
#set enum(spacing: 1.1em)   // room between list items so tall inline fractions don't collide
#set list(spacing: 1.1em)
#set table(stroke: 0.5pt + rgb("#c7ccd3"))

#let BLUE = rgb("#183B6F")     // brand navy — section bars + example/solution labels (white on navy 11.09:1)
#let ORANGE = rgb("#A94509")   // brand primary-700 — AA-safe deep orange for TEXT (5.93:1 on white; raw brand #F37021 is 2.94:1 and must never carry text)
#let RED = rgb("#DC2626")      // brand error-600
#let GREEN = rgb("#059669")    // brand success-600 (decoration only; small green text uses green-text #007942)

#show heading.where(level: 1): it => block(width: 100%, above: 0pt, below: 16pt,
  fill: gradient.linear(BLUE, rgb("#2C5AA0")), inset: (x: 14pt, y: 12pt), radius: 3pt,
  text(fill: white, weight: "bold", size: 19pt, it.body))
#show heading.where(level: 2): it => block(width: 100%, above: 18pt, below: 10pt,
  fill: BLUE, inset: (x: 10pt, y: 6pt), radius: 2pt,
  text(fill: white, weight: "bold", size: 12pt, it.body))
#show heading.where(level: 3): it => text(fill: ORANGE, weight: "bold", size: 12.5pt, it.body)
#show heading.where(level: 4): it => text(fill: BLUE, weight: "bold", size: 10.5pt, it.body)

#let examplebox(label, title, body) = block(width: 100%, breakable: true,
  fill: rgb("#EFF1F5"), stroke: 0.5pt + rgb("#CFDDF0"), radius: 4pt, inset: 10pt,
  above: 12pt, below: 12pt)[
  #block(below: 6pt)[#box(fill: BLUE, inset: (x: 6pt, y: 2pt), radius: 2pt,
    text(fill: white, weight: "bold", size: 8.5pt, label)) #h(0.4em) #strong[#title]]
  #body]

// rail = decorative left rule (raw brand token); labelcolor = AA-safe label text shade
#let notebox(label, rail, labelcolor, tint, body) = block(width: 100%, breakable: true, fill: tint,
    stroke: (left: 3pt + rail), inset: (left: 10pt, rest: 8pt), radius: (right: 4pt),
    above: 11pt, below: 11pt)[
    #text(fill: labelcolor, weight: "bold", size: 7.5pt, tracking: 0.5pt)[#upper(label)]
    #linebreak() #body]

#let solutionbox(body) = block(above: 4pt, below: 8pt)[
  #text(fill: BLUE, weight: "bold", size: 8.5pt)[Solution] #linebreak() #body]

#let figph(msg) = block(width: 100%, height: 60pt, fill: rgb("#f6f7f9"),
  stroke: (paint: rgb("#c7ccd3"), dash: "dashed"), radius: 4pt, inset: 10pt)[
  #align(center + horizon, text(fill: rgb("#889"), style: "italic", size: 9pt, msg))]

// Standardize inlined figure sizes: measure the natural CeTZ canvas, then scale to a
// consistent envelope (aspect-aware; see build_typst.py FIG_* constants). Unlike the
// print preamble, dimensions are FLOORED: in an editor a user can trim a figure to a
// degenerate 1-D shape (a bare line), and w/h or tw/w would then divide by zero.
#let _STD_W = 3.5
#let _WIDE_W = 5.6
#let _MAX_H = 3.4
#let _ASPECT_WIDE = 2.2
#let _UPSCALE_MAX = 1.15
#let stdfig(body) = context {
  let m = measure(body)
  let w = calc.max(m.width / 1in, 0.01)
  let h = calc.max(m.height / 1in, 0.01)
  let tw = if w / h > _ASPECT_WIDE { _WIDE_W } else { _STD_W }
  let s = calc.min(tw / w, _MAX_H / h, _UPSCALE_MAX)
  align(center, box(scale(x: s * 100%, y: s * 100%, reflow: true, body)))
}

#show figure: set block(breakable: false)
#set figure(gap: 8pt)
#show figure.caption: set text(size: 8.5pt, fill: rgb("#555"))


  
== 6.2#h(0.6em)Classification Using Machine Learning


  
    
=== Learning Outcomes


    By the end of this section, you should be able to:


    - Perform logistic regression on datasets and interpret the results.
- Perform k-means clustering on datasets.
- Define the concept of density-based clustering and use DBScan on datasets.
- Interpret the confusion matrix in clustering or classifying data.


  
  Classification problems come in many flavors. For example, suppose you are tasked with creating an algorithm that diagnoses heart disease. The input features may include biological sex, age, blood pressure data, and cholesterol levels. The output would be either: yes (diagnose heart disease) or no (do not diagnose heart disease). That is, your algorithm should classify patients as “yes” or “no” based on an array of features, or #emph[symptoms] in medical terminology. Logistic regression is one tool for classification when there are only two possible outputs. This is often called a #strong[binary (binomial) classification] problem. If the goal is to classify data into more than two classes or categories, then the problem is referred to as #strong[multiclass (multinomial) classification].


  Other kinds of classification problems involve finding two or more clusters in the data. A #strong[cluster] is a collection of data points that are closer to one another than to other data points, according to some definition of #emph[closeness]. For example, certain characteristics of music, such as tempo, instrumentation, overall length, and types of chord patterns used, can be used as features to classify the music into various genres, like rock, pop, country, and hip-hop. A supervised machine learning algorithm such as a decision tree (see Decision Trees) or random forest (see Other Machine Learning Techniques) may be trained and used to classify music into the various predefined genres. Alternatively, an unsupervised model such as k-means clustering could be used to group similar-sounding songs together without necessarily adhering to a concept of #emph[genre.]


  This section will focus on logistic regression techniques and clustering algorithms such as k-means and DBScan. We’ll use measures of accuracy including the confusion matrix.


  
    
=== Logistic Regression


    A #strong[logistic regression] model #math.equation(block: false, alt: "L open parenthesis X close parenthesis")[$L ( X )$] takes input vectors #math.equation(block: false, alt: "X")[$X$] (features) and produces an output of “yes” or “no.” Much like linear regression, the logistic regression model first fits a continuous function based on known, labeled data. However, instead of fitting a line to the data, a sigmoid function is used. A #strong[sigmoid function] is a specific type of function that maps any real-valued number to a value between 0 and 1. The term #emph[sigmoid] comes from the fact that the graph has a characteristic “S” shape (as shown in ), and #emph[sigma] #math.equation(block: false, alt: "open parenthesis σ close parenthesis")[$( σ )$] is the Greek letter that corresponds to our letter #emph[S]. We often use the notation  #math.equation(block: false, alt: "σ open parenthesis x close parenthesis")[$σ ( x )$] to denote a sigmoid function. The basic formula for the sigmoid is:


    #math.equation(block: true, alt: "σ open parenthesis x close parenthesis equals the fraction 1 over 1 plus e to the power minus x")[$σ ( x ) = frac(1, 1 + e^(− x))$]
    #figure(figph[A line graph.” The X axis ranges from -1 to 1 and a y axis in the center of the graph ranges from 0 to 1. The graph displays a sigmoid function that forms an S-curve that starts at (-1, 0), crosses the center Y axis at (0, 0.5), and ends at (1, 1).], alt: "A line graph.” The X axis ranges from -1 to 1 and a y axis in the center of the graph ranges from 0 to 1. The graph displays a sigmoid function that forms an S-curve that starts at (-1, 0), crosses the center Y axis at (0, 0.5), and ends at (1, 1).", caption: [Graph of a Typical Sigmoid Function])


    One advantage of the sigmoid over a linear function is that the values of the sigmoid stay close to 0 when #math.equation(block: false, alt: "x")[$x$] is negative and close to 1 when #math.equation(block: false, alt: "x")[$x$] is positive. This makes it ideal for classifying arbitrary inputs #math.equation(block: false, alt: "x")[$x$] into one of two classes, “yes” (1) or “no” (0). In fact, we can interpret the #emph[y]-axis as probability that an event occurs. The sigmoid function can be shifted right or left and compressed or expanded horizontally as necessary to fit the data. In the following formula, the parameters #math.equation(block: false, alt: "a")[$a$] and #math.equation(block: false, alt: "b")[$b$] control the position and shape of the sigmoid function, respectively.


    #math.equation(block: true, alt: "σ open parenthesis a plus b x close parenthesis equals the fraction 1 over 1 plus e to the power minus open parenthesis a plus b x close parenthesis")[$σ ( a + b x ) = frac(1, 1 + e^(− ( a + b x )))$]
    As complicated as this formula may seem at first, it is built up from simple elements. The presence of #math.equation(block: false, alt: "a plus b x")[$a + b x$]  suggests that there may be a linear regression #math.equation(block: false, alt: "y equals a plus b x")[$y = a + b x$] lurking in the background. Indeed, if we solve to find the inverse function of #math.equation(block: false, alt: "σ open parenthesis x close parenthesis")[$σ ( x )$], we can use it to isolate #math.equation(block: false, alt: "a plus b x")[$a + b x$], which in turns allows us to use standard linear regression techniques to build the logistic regression. Recall that the natural logarithm function, ln #math.equation(block: false, alt: "n")[$n$], is the inverse function for the exponential, #math.equation(block: false, alt: "e to the power x")[$e^(x)$]. In what follows, let #math.equation(block: false, alt: "p")[$p$] stand for   #math.equation(block: false, alt: "σ open parenthesis x close parenthesis")[$σ ( x )$].


    #math.equation(block: true, alt: "p equals σ open parenthesis x close parenthesis equals the fraction 1 over 1 plus e to the power minus x")[$p = σ ( x ) = frac(1, 1 + e^(− x))$]
    #math.equation(block: true, alt: "the fraction 1 over p, equals, 1 plus e to the power minus x; the fraction 1 over p minus 1, equals, e to the power minus x; e to the power minus x, equals, the fraction 1 over p minus 1 equals the fraction 1 minus p over p; minus x, equals, ln open parenthesis the fraction 1 minus p over p close parenthesis; x, equals, minus ln open parenthesis the fraction 1 minus p over p close parenthesis equals ln [ open parenthesis the fraction 1 minus p over p close parenthesis to the power minus 1 ] equals ln open parenthesis the fraction p over 1 minus p close parenthesis")[$frac(1, p) & = & 1 + e^(− x) \ frac(1, p) − 1 & = & e^(− x) \ e^(− x) & = & frac(1, p) − 1 = frac(1 − p, p) \ − x & = & "ln" ( frac(1 − p, p) ) \ x & = & − "ln" ( frac(1 − p, p) ) = "ln" [ attach(( frac(1 − p, p) ), t: − 1) ] = "ln" ( frac(p, 1 − p) )$]
    In the final line, we used the power property of logarithms, #math.equation(block: false, alt: "ln A to the power n equals n ln A")[$"ln " A^(n) = n " ln " A$], with power #math.equation(block: false, alt: "n equals minus 1")[$n = − 1$]. The function we’ve just obtained is called the logit function.


    #math.equation(block: true, alt: "logit open parenthesis p close parenthesis equals ln open parenthesis the fraction p over 1 minus p close parenthesis")[$"logit" ( p ) = "ln" ( frac(p, 1 − p) )$]
    The key property of the #strong[logit function] is that it is the inverse function of #math.equation(block: false, alt: "σ open parenthesis x close parenthesis")[$σ ( x )$]. Recall from algebra, if #math.equation(block: false, alt: "g open parenthesis x close parenthesis")[$g ( x )$] is the inverse function of #math.equation(block: false, alt: "f open parenthesis x close parenthesis")[$f ( x )$], then we have #math.equation(block: false, alt: "g open parenthesis f open parenthesis x close parenthesis close parenthesis equals x")[$g ( f ( x ) ) = x$] for all #math.equation(block: false, alt: "x")[$x$]. In our case, that means #math.equation(block: false, alt: "logit open parenthesis σ open parenthesis x close parenthesis close parenthesis equals x")[$"logit" ( σ ( x ) ) = x$], for all #math.equation(block: false, alt: "x")[$x$]. Now applying the logit function to the general sigmoid function, we can write:


    #math.equation(block: true, alt: "logit open parenthesis σ open parenthesis a plus b x close parenthesis close parenthesis equals a plus b x")[$"logit" #h(0.2em) ( σ ( a + b x ) ) = a + b x$]
    Thus, the logit function #emph[linearizes] the sigmoid. Let’s dig a little deeper into the logit function and what it measures.


    
      
==== Log-Odds and Discrete Logistic Regression


      Logistic regression is a prediction model based on the idea of odds. The #strong[#strong[odds]] that an event #emph[E] occurs is the ratio of the probability of #emph[E] occurring over the probability of #emph[E] not occurring. Let #math.equation(block: false, alt: "p equals P open parenthesis E close parenthesis")[$p = P ( E )$].


      #math.equation(block: true, alt: "Odds open parenthesis E close parenthesis equals the fraction P open parenthesis E close parenthesis over P open parenthesis E ′ close parenthesis equals the fraction P open parenthesis E close parenthesis over 1 minus P open parenthesis E close parenthesis equals the fraction p over 1 minus p")[$"Odds" ( E ) = frac(P ( E ), P ( E ′ )) = frac(P ( E ), 1 − P ( E )) = frac(p, 1 − p)$]
      Thus, the logit function is simply the logarithm of the odds of an event with a given probability #math.equation(block: false, alt: "p")[$p$]. This is why logit is often called #emph[log-odds]. While probabilities #math.equation(block: false, alt: "open parenthesis p close parenthesis")[$( p )$] range from 0 to 1, log-odds will range from #math.equation(block: false, alt: "minus infinity")[$− ∞$] (minus infinity) to #math.equation(block: false, alt: "infinity")[$∞$] (infinity).


      Suppose we would like to predict the likelihood of an event #math.equation(block: false, alt: "A")[$A$], given that another event #math.equation(block: false, alt: "B")[$B$] either occurs or does not occur. We will create a logistic regression model, #math.equation(block: false, alt: "f open parenthesis x close parenthesis equals the fraction 1 over 1 plus e to the power minus open parenthesis a plus b x close parenthesis")[$f ( x ) = frac(1, 1 + e^(− ( a + b x )))$], such that:


      - #math.equation(block: false, alt: "f open parenthesis 0 close parenthesis")[$f ( 0 )$] is the likelihood that #math.equation(block: false, alt: "A")[$A$] occurs when #math.equation(block: false, alt: "B")[$B$] does not occur, and
- #math.equation(block: false, alt: "f open parenthesis 1 close parenthesis")[$f ( 1 )$] is the likelihood that #math.equation(block: false, alt: "A")[$A$] occurs when #math.equation(block: false, alt: "B")[$B$] does occur.


      We obtain #math.equation(block: false, alt: "a")[$a$] and #math.equation(block: false, alt: "b")[$b$] using the logit function. If #math.equation(block: false, alt: "p sub 1")[$p_(1)$] is the probability that #math.equation(block: false, alt: "A")[$A$] occurs given that #math.equation(block: false, alt: "B")[$B$] does not occur (i.e., #math.equation(block: false, alt: "p sub 1 equals P open parenthesis A | B ′ close parenthesis")[$p_(1) = P ( A | B ′ )$]), using the notation of conditional probability (see Normal Continuous Probability Distributions), and if #math.equation(block: false, alt: "p sub 2")[$p_(2)$] is the probability that #math.equation(block: false, alt: "A")[$A$] occurs given that #math.equation(block: false, alt: "B")[$B$] does occur (i.e., #math.equation(block: false, alt: "p sub 2 equals P open parenthesis A | B close parenthesis")[$p_(2) = P ( A | B )$]), then


      #math.equation(block: false, alt: "a equals logit open parenthesis p sub 1 close parenthesis")[$a = "logit" #h(0.2em) ( p_(1) )$], and   #math.equation(block: false, alt: "b equals logit open parenthesis p sub 2 close parenthesis minus logit open parenthesis p sub 1 close parenthesis")[$b = "logit" #h(0.2em) ( p_(2) ) − "logit" #h(0.2em) ( p_(1) )$]


      This model is called a discrete logistic regression because the feature variable is either 1 or 0 (whether or not event #math.equation(block: false, alt: "B")[$B$] occurred).


      #examplebox("Example 1")[][
At a particular college, it has been found that in-state students are more likely to graduate than out-of-state students. Roughly 75% of in-state students end up graduating, while only about 30% of out-of-state students do so. Build a discrete logistic regression model based on this data.


        
        #solutionbox[
Let #math.equation(block: false, alt: "p sub 1 equals 0.3")[$p_(1) = 0.3$], which is the probability of graduating if the variable “in-state” is false, or 0. Let #math.equation(block: false, alt: "p sub 2 equals 0.75")[$p_(2) = 0.75$], which is the probability of graduating if the variable “in-state” is true, or 1.


          #math.equation(block: true, alt: "logit open parenthesis p sub 1 close parenthesis, equals, ln open parenthesis the fraction 0.3 over 1 minus 0.3 close parenthesis equals ln 0.4286 equals minus 0.847; logit open parenthesis p sub 2 close parenthesis, equals, ln open parenthesis the fraction 0.75 over 1 minus 0. 75 close parenthesis equals ln 3.000 equals 1.099")[$"logit" #h(0.2em) ( p_(1) ) & = & "ln " #h(0.2em) ( frac(0.3, 1 − 0.3) ) = "ln " 0.4286 = − 0.847 \ "logit" #h(0.2em) ( p_(2) ) & = & "ln " #h(0.2em) ( frac(0.75, 1 − 0. 75) ) = "ln " 3.000 = 1.099$]
          So #math.equation(block: false, alt: "a equals logit open parenthesis p sub 1 close parenthesis equals minus 0.847")[$a = "logit" #h(0.2em) ( p_(1) ) = − 0.847$], and #math.equation(block: false, alt: "b equals logit open parenthesis p sub 2 close parenthesis minus logit open parenthesis p sub 1 close parenthesis equals 1.099 minus open parenthesis minus 0.847 close parenthesis equals 1.946")[$b = "logit" #h(0.2em) ( p_(2) ) − "logit" #h(0.2em) ( p_(1) ) = 1.099 − ( − 0.847 ) = 1.946$]. Therefore, the model is:


          #math.equation(block: true, alt: "f open parenthesis x close parenthesis equals the fraction 1 over 1 plus e to the power minus open parenthesis minus 0.847 plus 1.946 x close parenthesis equals the fraction 1 over 1 plus e to the power 0.847 minus 1.946 x")[$f ( x ) = frac(1, 1 + e^(− ( − 0.847 + 1.946 x ))) = frac(1, 1 + e^(0.847 − 1.946 x))$]
          You can verify that #math.equation(block: false, alt: "f open parenthesis 0 close parenthesis equals 0.3")[$f ( 0 ) = 0.3$] and #math.equation(block: false, alt: "f open parenthesis 1 close parenthesis equals 0.75")[$f ( 1 ) = 0.75$].
]
]


      The previous example used only a single feature variable, “in-state,” and so the usefulness of the model is very limited. We shall see how to expand the model to include multiple inputs or features in Multiple Regression Techniques.


    
    
      
==== Maximum Likelihood and Continuous Logistic Regression


      The method of logistic regression can also be used to predict a yes/no response based on continuous input. Building the model requires finding values for the parameters of the sigmoid function that produce the most accurate results. To explain how this works, let’s first talk about likelihood. The #strong[likelihood] for a model #math.equation(block: false, alt: "f open parenthesis x close parenthesis")[$f ( x )$] is computed as follows. For each data point #math.equation(block: false, alt: "x sub k")[$x_(k)$], if the label of #math.equation(block: false, alt: "x sub k")[$x_(k)$] is 1 (or yes), then calculate #math.equation(block: false, alt: "f open parenthesis x sub k close parenthesis")[$f ( x_(k) )$]; if the label is 0 (or no), then calculate #math.equation(block: false, alt: "1 minus f open parenthesis x sub k close parenthesis")[$attach(1 − f ( x, b: k) )$]. The product of these individual results is the likelihood score for the model. We want to find the model that has the maximum likelihood.


      #examplebox("Example 2")[][
Consider the models


          #math.equation(block: true, alt: "f open parenthesis x close parenthesis equals the fraction 1 over 1 plus e to the power 1.6 minus 1.1 x , g open parenthesis x close parenthesis equals the fraction 1 over 1 plus e to the power 2.1 minus 0.8 x")[$f ( x ) = frac(1, 1 + e^(1.6 − 1.1 x)) , g ( x ) = frac(1, 1 + e^(2.1 − 0.8 x))$]
          Suppose there are four data points, #math.equation(block: false, alt: "A equals 1")[$A = 1$], #math.equation(block: false, alt: "B equals 2")[$B = 2$], #math.equation(block: false, alt: "C equals 3")[$C = 3$], #math.equation(block: false, alt: "D equals 4")[$D = 4$], and both #math.equation(block: false, alt: "A")[$A$] and #math.equation(block: false, alt: "B")[$B$] are known to have label 0, while #math.equation(block: false, alt: "C")[$C$] and #math.equation(block: false, alt: "D")[$D$] have label 1. Find the likelihood scores of each model and determine which model is a better fit for the given data.


        
        #solutionbox[
For model #math.equation(block: false, alt: "f")[$f$], the predicted likelihoods are   #math.equation(block: false, alt: "f open parenthesis 1 close parenthesis equals 0.378")[$f ( 1 ) = 0.378$], #math.equation(block: false, alt: "f open parenthesis 2 close parenthesis equals 0.646")[$f ( 2 ) = 0.646$],   #math.equation(block: false, alt: "f open parenthesis 3 close parenthesis equals 0.846")[$f ( 3 ) = 0.846$], and #math.equation(block: false, alt: "f open parenthesis 4 close parenthesis equals 0.943")[$f ( 4 ) = 0.943$]. Since #math.equation(block: false, alt: "x equals 1")[$x = 1$] and 2 have label 0, while #math.equation(block: false, alt: "x equals 3")[$x = 3$] and 4 have label 1, we compute:


          #math.equation(block: true, alt: "open parenthesis 1 minus 0.378 close parenthesis open parenthesis 1 minus 0.646 close parenthesis open parenthesis 0.846 close parenthesis open parenthesis 0.943 close parenthesis equals 0.176")[$( 1 − 0.378 ) ( 1 − 0.646 ) ( 0.846 ) ( 0.943 ) = 0.176$]
          For model #math.equation(block: false, alt: "g")[$g$], the predicted likelihoods are #math.equation(block: false, alt: "g open parenthesis 1 close parenthesis equals 0.214")[$g ( 1 ) = 0.214$], #math.equation(block: false, alt: "g open parenthesis 2 close parenthesis equals 0.378")[$g ( 2 ) = 0.378$], #math.equation(block: false, alt: "g open parenthesis 3 close parenthesis equals 0.574")[$g ( 3 ) = 0.574$], and #math.equation(block: false, alt: "g open parenthesis 4 close parenthesis equals 0.75")[$g ( 4 ) = 0.75$].


          #math.equation(block: true, alt: "open parenthesis 1 minus 0.214 close parenthesis open parenthesis 1 minus 0.378 close parenthesis open parenthesis 0.574 close parenthesis open parenthesis 0.75 close parenthesis equals 0.211")[$( 1 − 0.214 ) ( 1 − 0.378 ) ( 0.574 ) ( 0.75 ) = 0.211$]
          Thus, model #math.equation(block: false, alt: "g")[$g$] is the better of the two models on this data.
]
]


      How does one find the model that has maximum likelihood? The exact method for finding the coefficients #math.equation(block: false, alt: "a")[$a$] and #math.equation(block: false, alt: "b")[$b$] for the model #math.equation(block: false, alt: "σ open parenthesis a plus b x close parenthesis")[$σ ( a + b x )$] falls outside the scope of this text. Fortunately, there are software packages available that can be used to do the work for you.


      For example, the same university ran a study on the effect of first-year GPA on completion of a college degree. A small sample of the data is shown. What is the model that best fits this data?


      #figure(table(
  columns: 2, align: left, inset: 6pt,
  table.header([#strong[GPA]], [#strong[Completed College?]]),
  [1.5], [N],
  [2.4], [Y],
  [3.4], [Y],
  [2.1], [N],
  [2.5], [Y],
  [0.8], [N],
  [2.9], [N],
  [4.0], [Y],
  [2.3], [Y],
  [2.1], [N],
  [3.6], [Y],
  [0.5], [N],
))


      shows the data in graphical form.


      #figure(figph[A line graph with an X axis that ranges from 0 to 4 and a Y axis that ranges from 0 to 1. There are five data points at 0 on the Y axis (0, 0.5,) (0, 0.8), (0, 1.5), (0, 2.1), (0, 2.9) and six data points at 1 on the Y axis at (2.3, 1), (2.4, 1), ( 2.5, 1), (3.4, 1), (3.6, 1), ( 4.0, 1).], alt: "A line graph with an X axis that ranges from 0 to 4 and a Y axis that ranges from 0 to 1. There are five data points at 0 on the Y axis (0, 0.5,) (0, 0.8), (0, 1.5), (0, 2.1), (0, 2.9) and six data points at 1 on the Y axis at (2.3, 1), (2.4, 1), ( 2.5, 1), (3.4, 1), (3.6, 1), ( 4.0, 1).", caption: [Graph of GPA (#emph[x]-axis) and College Completion (#emph[y]-axis). This uses values 0 for non-completion and 1 for completion. The graph suggests that higher GPAs tend to predict completion of a degree.])


      The logistics model (found using computer software) that maximizes the likelihood is:


      #math.equation(block: true, alt: "the fraction 1 over 1 plus e to the power 2.826 minus 1.197 x")[$frac(1, 1 + e^(2.826 − 1.197 x))$]
      The logistics regression model that best fits the data is shown in  along with the data points colored green for predicted completion and red for predicted non-completion.


      #figure(figph[A line graph with an X axis that ranges from 0 to 4 and a Y axis that ranges from 0 to 1. There are five red data points (0, 0.5,) (0, 0.8), (0, 1.5), (0, 2.1), and (2.3, 1) and six green data points (0, 2.9), (2.4, 1), ( 2.5, 1), (3.4, 1), (3.6, 1), ( 4.0, 1). Data points are colored green for predicted completion and red for predicted non-completion.], alt: "A line graph with an X axis that ranges from 0 to 4 and a Y axis that ranges from 0 to 1. There are five red data points (0, 0.5,) (0, 0.8), (0, 1.5), (0, 2.1), and (2.3, 1) and six green data points (0, 2.9), (2.4, 1), ( 2.5, 1), (3.4, 1), (3.6, 1), ( 4.0, 1). Data points are colored green for predicted completion and red for predicted non-completion.", caption: [Logistic Regression Model Fitting the Data in])


    
    
      
==== Logistic Regression in Python


      A sample of Python code for logistic regression, using the data from , appears next. We use the #strong[LogisticRegression] function found in the Python library #link("https://openstax.org/r/scikit")[#strong[sklearn.linear\_model]].


      #notebox("Note", rgb("#8a94a6"), rgb("#556666"), rgb("#f7f8fa"))[
\# Import libraries
      from sklearn.linear\_model import LogisticRegression
      
      \# Define input and output data
      data = \[\[1.5, 0\], \[2.4, 1\], \[3.4, 1 \], \[2.1, 0 \], \[2.5, 1 \], 
       \[0.8, 0 \], \[2.9, 0 \], \[4.0, 1\], \[2.3, 1\], \[2.1, 0\], \[3.6, 1 \], \[0.5, 0\]\]
      
      \# Separate features (X) and labels (y)
      X = \[\[row\[0\]\] for row in data\]  \# Feature
      y = \[row\[1\] for row in data\]    \# Label
      
      \# Build the logistic model
      Lmodel = LogisticRegression()
      Lmodel.fit(X,y)
      
      \# Display the coefficients of the logistic regression
      (a,b) = (Lmodel.intercept\_\[0\], Lmodel.coef\_\[0,0\])
      print("Coefficients: ", (round(a,3),round(b,3)))
      
      \# Display the accuracy of the model
      s = Lmodel.score(X,y)
      print("Accuracy: ", round(s, 3))
    
        The resulting output will look like this:


        Coefficients: (-2.826, 1.197) #linebreak() Accuracy: 0.833
]


      The code produces the coefficients that can be plugged in for #math.equation(block: false, alt: "a")[$a$] and #math.equation(block: false, alt: "b")[$b$] in the model:


      #math.equation(block: true, alt: "σ open parenthesis minus 2.826 plus 1.197 x close parenthesis equals the fraction 1 over 1 plus e to the power minus open parenthesis minus 2.826 plus 1.197 x close parenthesis equals the fraction 1 over 1 plus e to the power 2.826 minus 1.197 x")[$σ ( − 2.826 + 1.197 x ) = frac(1, 1 + e^(− ( − 2.826 + 1.197 x ))) = frac(1, 1 + e^(2.826 − 1.197 x))$]
    
  
  
    
=== k-Means Clustering


    Clustering algorithms perform like digital detectives, uncovering patterns and structure in data. The goal is to find and label groups of data points that are close to one another. This section develops methods for grouping data (clustering) that incorporate machine learning.


    The #strong[k-means clustering algorithm] can classify or group similar data points into clusters or categories without prior knowledge of what those categories might be (i.e., unsupervised learning). However, the user must provide a guess as to the value of #math.equation(block: false, alt: "k")[$k$], the number of clusters to expect. (Note: You may get around this issue by iterating the algorithm through many values of #math.equation(block: false, alt: "k")[$k$] and evaluating which #emph[k]-value produced the best results. The trade-off is the time it takes to run the algorithm using multiple values of #math.equation(block: false, alt: "k")[$k$] and interpreting the results using. The so-called elbow method is often used to find the best value of #math.equation(block: false, alt: "k")[$k$], but it is beyond the scope of this text to explain this method.) The user must also provide initial guesses as to where the centroids would be. We shall see by example how these choices may affect the algorithm.


    
      
==== k-Means Clustering Fundamentals


      The basic idea behind k-means is quite intuitive: Each cluster must be centered around its centroid in such a way that minimizes distance from the data in the cluster to its centroid. The trick is in finding those centroids! Briefly, a #strong[centroid] of a set of data points is defined as the arithmetic mean of the data, regarded as vectors.


      #math.equation(block: true, alt: "Centroid equals the fraction 1 over N ∑ i equals 1 N X sub i")[$"Centroid" #h(0.2em) = frac(1, N) limits(∑)^(N)_(i = 1) X_(i)$]
      The basic steps of the k-means clustering algorithm are as follows:


      + Guess the locations of #math.equation(block: false, alt: "k")[$k$] cluster centers.
+ Compute distances from each data point to each cluster center. Data points are assigned to the cluster whose center is closest to them.
+ Find the centroid of each cluster.
+ Repeat steps 2 and 3 until the centroids stabilize—that is, the new centroids are the same (or approximately the same, within some given tolerance) as the centroids from the previous iteration.


      If the clustering was done well, then clusters should be clearly separated. The #strong[silhouette score] is a measure of how well-separated the clusters are and is defined as follows:


      Let #math.equation(block: false, alt: "a open parenthesis i close parenthesis")[$a ( i )$] be the mean distance between point #math.equation(block: false, alt: "i")[$i$] and all other data points in the same cluster as point #math.equation(block: false, alt: "i")[$i$].


      For each cluster #math.equation(block: false, alt: "J")[$J$] not containing point #math.equation(block: false, alt: "i")[$i$], let #math.equation(block: false, alt: "b sub J open parenthesis i close parenthesis")[$b_(J) ( i )$] be the mean distance between point #math.equation(block: false, alt: "i")[$i$] and all data points in cluster #math.equation(block: false, alt: "J")[$J$]. Then, let #math.equation(block: false, alt: "b open parenthesis i close parenthesis")[$b ( i )$] be the minimum of all the values of #math.equation(block: false, alt: "b sub J open parenthesis i close parenthesis")[$b_(J) ( i )$].


      Then the silhouette score is found by the formula


      #math.equation(block: true, alt: "S equals the fraction 1 over N ∑ i equals 1 N the fraction b open parenthesis i close parenthesis minus a open parenthesis i close parenthesis over max { a open parenthesis i close parenthesis , b open parenthesis i close parenthesis }")[$S = frac(1, N) limits(∑)^(N)_(i = 1) frac(b ( i ) − a ( i ), "max" #h(0.2em) \{ a ( i ) "," b ( i ) \})$]
      Values of #math.equation(block: false, alt: "S")[$S$] lie between #math.equation(block: false, alt: "minus")[$−$]1 and 1, with values close to 1 indicating well-separated clusters, values close to 0 indicating clusters are ambiguous, and values close to #math.equation(block: false, alt: "minus")[$−$]1 indicating poor clustering with points assigned to clusters arbitrarily. Fortunately, statistical software packages that can do k-means clustering will be able to compute the silhouette score for you.


      #examplebox("Example 3")[][
An outbreak of fungus has affected a small garden area. The fungus originated underground, but it produced mushrooms that can easily be spotted. From a photo of the area, it is surmised that there are three clusters of mushrooms. Use k-means to classify the mushrooms into three groups.


          #figure(figph[A scatterplot with an X axis labeled “East” that ranges from 20 to 160 and a Y axis labeled “North” that ranges from 20 to 90. There are 19 data points that represent locations of mushrooms in a garden and represent three apparent clusters with a couple of outliers.], alt: "A scatterplot with an X axis labeled “East” that ranges from 20 to 160 and a Y axis labeled “North” that ranges from 20 to 90. There are 19 data points that represent locations of mushrooms in a garden and represent three apparent clusters with a couple of outliers.", caption: [The Locations of Mushrooms in a Garden. There are three apparent clusters.])


        
        #solutionbox[
The location data consists of 19 points indicating position in feet east and north of a fixed point.


          (101.7, 51.1), (109.2, 66.7), (33.2, 54.1), (39.6, 53.6), (66.8, 70.5), (119.9, 42.4),


          (26.9, 33.1), (55.5, 71.7), (75, 82.4), (108.6, 71.4), (62.1, 83.5), (119.8, 69.1),


          (25.9, 30.2), (37.5, 46.4), (24, 49.1), (75.5, 77.8), (101.5, 51), (110.5, 72),


          (39.5, 55.8)


          First, choose three potential centroids. To make things simple, we will just pick the first three data points (101.7, 51.1), (109.2, 66.7), and (33.2, 54.1), which we will color red, green, and blue, respectively. Now color all the data points red, green, and blue according to which of these three points it lies closest to. The initial grouping is not perfect, as you can see in .


          #figure(figph[A scatterplot with an X axis labeled “East” that ranges from 20 to 160 and a Y axis labeled “North” that ranges from 20 to 90. There are 19 data points that represent locations of mushrooms in a garden. The points are colored blue, green, or red according to whether the data points lie closest to one of three centroids (101.7, 51.1), (109.2, 66.7) or (33.2, 54.1), which are represent by larger dots.], alt: "A scatterplot with an X axis labeled “East” that ranges from 20 to 160 and a Y axis labeled “North” that ranges from 20 to 90. There are 19 data points that represent locations of mushrooms in a garden. The points are colored blue, green, or red according to whether the data points lie closest to one of three centroids (101.7, 51.1), (109.2, 66.7) or (33.2, 54.1), which are represent by larger dots.", caption: [Graph with Large Dots Indicating the Initial Choices of Centroids. Data points are colored according to which centroid they are closest to.])


          Next, we find the centroids of each of the colored groups. In practice, this is all done automatically by software, but we will work out one of the centroids explicitly here. Let’s work out the green centroid.


          Current green points: (75, 82.4), (75.5, 77.8), (108.6, 71.4), (109.2, 66.7), (110.5, 72), (119.8, 69.1).


          Green centroid #emph[x]-value: #math.equation(block: false, alt: "x minus equals the fraction 1 over 6 open parenthesis 75 plus 75.5 plus 108.6 plus 109.2 plus 110.5 plus 119.8 close parenthesis equals 99.8")[$limits(x)^(−) = frac(1, 6) ( 75 + 75.5 + 108.6 + 109.2 + 110.5 + 119.8 ) = 99.8$]


          Green centroid #emph[y]-value:   #math.equation(block: false, alt: "y minus equals the fraction 1 over 6 open parenthesis 82.4 plus 77.8 plus 71.4 plus 66.7 plus 72 plus 69.1 close parenthesis equals 73.2")[$limits(y)^(−) = frac(1, 6) ( 82.4 + 77.8 + 71.4 + 66.7 + 72 + 69.1 ) = 73.2$]


          Thus, the green centroid is located at (99.8, 73.2).


          Now plot all the centroids. You might imagine the previous potential centroids have moved to positions that are closer to the true centroids of the visible groups. The red point moves to (107.7, 48.2); the green point moves to (99.8, 73.2); and the blue point moves to (41.1, 54.8). Note that (aside from the initial choices) the centroids are not typically points in the dataset. Recolor the points of the dataset with respect to these new centroids, as shown.


          #figure(figph[A scatterplot with an X axis labeled “East” that ranges from 20 to 160 and a Y axis labeled “North” that ranges from 20 to 90. There are 19 data points that represent locations of mushrooms in a garden. The points are colored blue, green, or red according to whether the data points lie closest to one of three centroids (107.7, 48.2), (99.8, 73.2) or (41.1, 54.8), which are represent by larger dots.], alt: "A scatterplot with an X axis labeled “East” that ranges from 20 to 160 and a Y axis labeled “North” that ranges from 20 to 90. There are 19 data points that represent locations of mushrooms in a garden. The points are colored blue, green, or red according to whether the data points lie closest to one of three centroids (107.7, 48.2), (99.8, 73.2) or (41.1, 54.8), which are represent by larger dots.", caption: [Graph with New Centroids. After the second iteration of the k-means algorithm, the centroids have moved closer to the centers of the main clusters. However, no points changed location, so the algorithm stabilizes at this step.])


          Run the algorithm a number of times until the coloring stabilizes. In this case, the next iteration is stable. Are the clusters where you thought they might end up? This example highlights the fact that k-means may be quite sensitive to initial choice of centroids.


          Running the algorithm again with different initial centroids may produce different clusters, as  shows.


          #figure(figph[A scatterplot with an X axis labeled “East” that ranges from 20 to 160 and a Y axis labeled “North” that ranges from 20 to 90. There are 19 data points that represent locations of mushrooms in a garden. The points are colored blue, green, or red according to whether the data points lie closest to one of three centroids (110.2, 60.5), (67.0, 77.2), and (32.4, 46.0), which are represent by larger dots.], alt: "A scatterplot with an X axis labeled “East” that ranges from 20 to 160 and a Y axis labeled “North” that ranges from 20 to 90. There are 19 data points that represent locations of mushrooms in a garden. The points are colored blue, green, or red according to whether the data points lie closest to one of three centroids (110.2, 60.5), (67.0, 77.2), and (32.4, 46.0), which are represent by larger dots.", caption: [Graph with Different Initial Centroids. Using initial centroids of (110, 70), (60, 80), and (40, 40), the k-means algorithm produces a different clustering. The final position of the centroids are (110.2, 60.5), (67.0, 77.2), and (32.4, 46.0).])


          Clustering algorithms such as k-means also allow one to make predictions about new data. For example, suppose a new mushroom sprouted up at location (90, 50). To classify this point as red, green, or blue, simply find out which centroid is closest to the point.


          Recall that the distance between two points #math.equation(block: false, alt: "open parenthesis x sub 1 , y sub 1 close parenthesis")[$( x_(1) , y_(1) )$] and #math.equation(block: false, alt: "open parenthesis x sub 2 , y sub 2 close parenthesis")[$( x_(2) , y_(2) )$] is found using the formula


          #math.equation(block: true, alt: "d equals the square root of open parenthesis x sub 2 minus x sub 1 close parenthesis squared plus open parenthesis y sub 2 minus y sub 1 close parenthesis squared")[$d = sqrt(attach(( x_(2) − x_(1) ), t: 2) + attach(( y_(2) − y_(1) ), t: 2))$]
          Distance to red centroid: #math.equation(block: false, alt: "the square root of open parenthesis 90 minus 110.2 close parenthesis squared plus open parenthesis 50 minus 60.5 close parenthesis squared equals 22.77")[$sqrt(attach(( 90 − 110.2 ), t: 2) + attach(( 50 − 60.5 ), t: 2)) = 22.77$]


          Distance to green centroid: #math.equation(block: false, alt: "the square root of open parenthesis 90 minus 67 close parenthesis squared plus open parenthesis 50 minus 77.2 close parenthesis squared equals 35.62")[$sqrt(attach(( 90 − 67 ), t: 2) + attach(( 50 − 77.2 ), t: 2)) = 35.62$]


          Distance to blue: #math.equation(block: false, alt: "the square root of open parenthesis 90 minus 32.4 close parenthesis squared plus open parenthesis 50 minus 46 close parenthesis squared equals 57.74")[$sqrt(attach(( 90 − 32.4 ), t: 2) + attach(( 50 − 46 ), t: 2)) = 57.74$]


          Thus, the new mushroom should be included in the red group.
]
]


      Among clustering algorithms, k-means is rather simple, easy to implement, and very fast. However, it suffers from many drawbacks, including sensitivity to initial choice of centroids. Furthermore, k-means does detect clusters that are not simply “blobs.” In the next section, we will discuss another clustering algorithm that generally performs better on a wider range of data.


    
    
      
==== k-Means Clustering in Python


      Our previous example used a sample of only 19 data points selected from the 115 data points found in the file #link("https://openstax.org/r/FungusLocations")[#strong[FungusLocations.csv]]. Here is the Python code to produce a k-means clustering from a dataset:


      #notebox("Note", rgb("#8a94a6"), rgb("#556666"), rgb("#f7f8fa"))[
\# Import libraries
      import pandas as pd  \#\# for dataset management
      import matplotlib.pyplot as plt \#\# for data visualization
      from sklearn.cluster import KMeans
      from sklearn.metrics import silhouette\_score
      
      \# Read data
      data = pd.read\_csv('FungusLocations.csv').dropna()
      
      \# Build K-means model with 3 clusters
      km = KMeans(n\_clusters=3, n\_init='auto')
      km.fit(data)
      
      \# Calculate silhouette score
      silhouette\_avg = silhouette\_score(data, km.labels\_)
      print("Silhouette Score:", round(silhouette\_avg,2))
      
      \# Visualize the result of k-Means clustering using matplotlib
      plt.scatter(data\['East'\], data\['North'\], c=km.labels\_, cmap='viridis')
      plt.xlabel('East')
      plt.ylabel('North')
      plt.title('k-Means Clustering Result')
      plt.show()
      
        The resulting output will look like this:


        Silhouette Score: 0.77
        #figure(figph[A scatterplot labeled “k-Means Clustering Result” with an X axis labeled “East” that ranges from 20 to 160 and a Y axis labeled “North” that ranges from 20 to 90. There are 19 data points that represent locations of mushrooms in a garden. The points are shown in three clusters colored yellow, purple, and green.], alt: "A scatterplot labeled “k-Means Clustering Result” with an X axis labeled “East” that ranges from 20 to 160 and a Y axis labeled “North” that ranges from 20 to 90. There are 19 data points that represent locations of mushrooms in a garden. The points are shown in three clusters colored yellow, purple, and green.", caption: none)
]


      With a silhouette score of about 0.77, the separation of data points into three clusters seems appropriate.


      #notebox("Note", rgb("#8a94a6"), rgb("#556666"), rgb("#f7f8fa"))[
#emph[Centroids versus Other Means]


        Instead of using #emph[centroids] (which are computed using arithmetic means, or averages), other means may be employed, including medians, trimmed averages, geometric or harmonic means, and many others. These algorithms may have different names such as #emph[k-medians], #emph[k-medoids], etc. More information can be found at #link("https://openstax.org/r/neptune")[Neptune] and #link("https://openstax.org/r/tidymodels")[Tidymodels articles]..
]


    
  
  
    
=== Density-Based Clustering (DBScan)


    Suppose you need to classify the data in  into two clusters.


    #figure(figph[A density-based clustering scan with an X axis that ranges from -15 to 15 and a Y axis that ranges from -15 to 15. Blue dots form a larger outer circle cluster and a smaller inner cluster in the center of the circle.], alt: "A density-based clustering scan with an X axis that ranges from -15 to 15 and a Y axis that ranges from -15 to 15. Blue dots form a larger outer circle cluster and a smaller inner cluster in the center of the circle.", caption: [Dataset with an Inner and Outer Ring])


    How might you do this? There seems to be a clear distinction between points on the inner ring and points on the outer ring. While k-means clustering would not do well on this dataset, a #strong[density-based clustering algorithm] such as DBScan might perform better.


    
      
==== DBScan Fundamentals


      The #strong[DBScan algorithm] (DBScan is an acronym that stands for #emph[density-based spatial clustering of applications with noise]) works by locating core data points that are part of dense regions of the dataset and expanding from those cores by adding neighboring points according to certain closeness criteria.


      First, we define two parameters: #math.equation(block: false, alt: "K")[$K$], and #math.equation(block: false, alt: "r")[$r$]. A point in the dataset is called a core point if it has at least #math.equation(block: false, alt: "K")[$K$] neighbors, counting itself, within a distance of #math.equation(block: false, alt: "r")[$r$] from itself.


      The steps of DBScan are as follows:


      + Choose a starting point that is a core point and assign it to a cluster. (Note: If there are no core points in the dataset at all, then the parameters #math.equation(block: false, alt: "K")[$K$] and #math.equation(block: false, alt: "r")[$r$] will need to be adjusted.)
+ Add all core points that are close (less than a distance of #math.equation(block: false, alt: "r")[$r$]) to the first point into the first cluster. Keep adding to the cluster in this way until there are no more core points close enough to the first cluster.
+ If there are any core points not assigned to a cluster, choose one to start a new cluster and repeat step 2 until all core points are assigned to clusters.
+ For each non-core point, check to see if it is close to any core points. If so, then add it to the cluster of the closest core point. If a non-core point is not close enough to any core points, it will not be added to any cluster. Such a point is regarded as an outlier or noise.


      Since the values of #math.equation(block: false, alt: "K")[$K$] and #math.equation(block: false, alt: "r")[$r$] are so important to the performance of DBScan, the algorithm is typically run multiple times on the same dataset using a range of #math.equation(block: false, alt: "K")[$K$] values and #math.equation(block: false, alt: "r")[$r$] values. If there are a large number of outliers and/or very few clusters, then the model may be underfitting. Conversely, a very small number of outliers (or zero outliers) and/or many clusters present in the model may be a sign of overfitting.


      Let’s work through a small example using DBScan with #math.equation(block: false, alt: "K equals 3")[$K = 3$] and #math.equation(block: false, alt: "r equals 2.1")[$r = 2.1$] to classify the following dataset into clusters. (Points are labeled by letters for convenience.)


      A: (3.5, 5.8)	B: (4.5, 3.5)	C: (5.0, 5.0)	D: (6.0, 6.0)	E: (6.0, 3.0)	F: (6.9, 4.7)


      G: (9.0, 4.0)	H: (11.0, 5.0)	J: (12.0, 4.0)	K: (12.4, 5.4)	L: (13.2, 2.5)


      Note: In practice, DBScan would never be done by hand, but the steps shown as follows illustrate how the algorithm works.  displays the data.


      #figure(figph[A scatterplot with an X axis that ranges from 0 to 16 and a Y axis that ranges from 0 to 8. There are 11 data points labeled A through L randomly scattered.], alt: "A scatterplot with an X axis that ranges from 0 to 16 and a Y axis that ranges from 0 to 8. There are 11 data points labeled A through L randomly scattered.", caption: [Original Dataset with 11 Unlabeled Points.])


      The core points are those that have at least #math.equation(block: false, alt: "K")[$K$] neighbors (or #math.equation(block: false, alt: "K minus 1")[$K − 1$], not counting the point itself) within a distance of #math.equation(block: false, alt: "r")[$r$]. The easiest way to determine core points would be to draw circles of radius 2.1 centered at each point, as shown in , and then count the number of points within each circle.


      #figure(figph[A scatterplot with an X axis that ranges from 0 to 16 and a Y axis that ranges from 0 to 8. There are 11 data points labeled A through L randomly scattered. There are 11 blue circles, each centered around a data point. The circles overlap and intersect illustrating which points are core or not core.], alt: "A scatterplot with an X axis that ranges from 0 to 16 and a Y axis that ranges from 0 to 8. There are 11 data points labeled A through L randomly scattered. There are 11 blue circles, each centered around a data point. The circles overlap and intersect illustrating which points are core or not core.", caption: [Circles of Radius Centered at Each Point. Drawing circles of radius #math.equation(block: false, alt: "r")[$r$] around each data point, it is easier to see which points are core and which are not core.])


      Based on , point #emph[A] has only two neighbors, #emph[A] and #emph[C], so #emph[A] is not a core point. On the other hand, point #emph[B] has three neighbors, #emph[B], #emph[C], and #emph[E], so #emph[B] is core. The full list of core points is #emph[B], #emph[C], #emph[D], #emph[E], #emph[F], #emph[H], #emph[J], and #emph[K]. So #emph[B] may be chosen to start the first cluster. Then, building from #emph[B], points #emph[C] and #emph[E] are added next. Next, #emph[D], #emph[E], and #emph[F] are added to the growing first cluster. Neither #emph[A] nor #emph[G] are added at this stage since they are not core points. Having found all core points that are in the first cluster, let #emph[H] start a second cluster. The only core points left are close enough to #emph[H] that the second cluster is #emph[H], #emph[J], and #emph[K].


      Finally, we try to place the non-core points into clusters. Point #emph[A] is close enough to #emph[C] to make it into the first cluster. Point #emph[L] is close enough to #emph[J] to become a part of the second cluster. Point #emph[G] does not have any core points near enough, so #emph[G] will be un-clustered and considered noise. Thus, cluster 1 consists of #emph[{A, B, C, D, E, F}], and cluster 2 consists of #emph[{H, J, K, L}]. The results of DBScan clustering are shown.


      #figure(figph[A scatterplot with an X axis that ranges from 0 to 16 and a Y axis that ranges from 0 to 8. There are 11 data points labeled A through L. Data points A-F are blue, G is black, and H-L are red, representing a DBScan of two clusters. Cluster 1 is in blue and Cluster 2 is in red.], alt: "A scatterplot with an X axis that ranges from 0 to 16 and a Y axis that ranges from 0 to 8. There are 11 data points labeled A through L. Data points A-F are blue, G is black, and H-L are red, representing a DBScan of two clusters. Cluster 1 is in blue and Cluster 2 is in red.", caption: [Graph Showing Results of DBScan Clustering. Cluster 1 is in blue and Cluster 2 is in red.])


    
    
      
==== DBScan in Python


      The Python library #link("https://openstax.org/r/stable")[#strong[sklearn.cluster]] has a module named #strong[DBSCAN]. Here is how it works on the dataset #link("https://openstax.org/r/DBScanExample")[#strong[DBScanExample.csv]].


      #notebox("Note", rgb("#8a94a6"), rgb("#556666"), rgb("#f7f8fa"))[
\# Import libraries
      import pandas as pd \#\# for dataset management
      from sklearn.cluster import DBSCAN
      
      \# Read data
      data = pd.read\_csv('DBScanExample.csv').dropna()
      
       \# Run DBScan
      db = DBSCAN(eps=1.3, min\_samples=5).fit(data)
      labels = db.labels\_
      
      \# Find the number of clusters and points of noise
      n\_clusters = len(set(labels)) - (1 if -1 in labels else 0)
      n\_noise = list(labels).count(-1)
      print("Number of clusters: ", n\_clusters)
      print("Points of noise:", n\_noise)
    
        The resulting output will look like this:


        Number of clusters: 2 #linebreak() Points of noise: 8
]


      #notebox("Note", rgb("#8a94a6"), rgb("#556666"), rgb("#f7f8fa"))[
\# Visualize clusters
      import matplotlib.pyplot as plt \#\# for data visualization
      
      plt.scatter(data\['X'\], data\['Y'\], c=db.labels\_)
    
        The resulting output will look like this:


        #figure(figph[A density-based clustering scan with an X axis that ranges from -10 to 10 and a Y axis that ranges from -10 to 10. Yellow dots form a larger outer circle cluster and green dots form a smaller inner cluster in the center of the circle. Eight purple dots appear as outliers inside and outside the circle.], alt: "A density-based clustering scan with an X axis that ranges from -10 to 10 and a Y axis that ranges from -10 to 10. Yellow dots form a larger outer circle cluster and green dots form a smaller inner cluster in the center of the circle. Eight purple dots appear as outliers inside and outside the circle.", caption: none)
]


    
  
  
    
=== The Confusion Matrix


    When training and testing any algorithm that performs classification, such as logistic regression models, k-means clustering, or DBScan, it is important to measure accuracy and identify error. The #strong[confusion matrix] is one way of quantifying and visualizing the effectiveness of a classification model. This concept will be illustrated by the following example.


    Suppose you have trained a model that classifies images of plants into one of three categories: flowers, trees, and grasses. When you test the model on 100 additional images of plants, you discover that it correctly identified flowers, trees, and grasses the majority of the time, but there were also quite a few mistakes.  displays the results in a confusion matrix:


    #figure(table(
  columns: 4, align: left, inset: 6pt,
  table.header([], [#strong[Identified as a Flower]], [#strong[Identified as a Tree]], [#strong[Identified as a Grass]]),
  [Is a flower], [23], [3], [9],
  [Is a tree], [2], [32], [0],
  [Is a grass], [12], [1], [18],
))


    Thus, we can see that 23 flowers, 32 trees, and 18 grasses were identified correctly, giving the model an accuracy of #math.equation(block: false, alt: "the fraction 23 plus 32 plus 18 over 100 equals 0.73")[$frac(23 + 32 + 18, 100) = 0.73$], or 73%. Note: The terms TP, TN, FP, and FN introduced in What Is Machine Learning? and defined again in the next example do not apply when classifying data into three or more categories; however, the confusion matrix gives more detailed information about the mistakes. Where there are relatively low numbers off the main diagonal, the model performed well. Conversely, higher numbers off the main diagonal indicated greater rates of misidentification, or confusion. In this example, when presented with a picture of a flower, the model did a pretty good job of identifying it as a flower but would mistake flowers for grasses at a higher rate than mistaking them for trees. The model did very well identifying trees as trees, never mistaking them for grasses and only twice mistaking them for flowers. On the other hand, the model did not do so well identifying grasses. Almost 40% of the images of grasses were mistaken for flowers! With this information in hand, you could go back to your model and adjust parameters in a way that may address these specific issues. There may be a way to tweak the model that helps with identifying grasses in particular. Or you may decide to take extra steps whenever an image is identified as a flower, perhaps running the image through additional analysis to be sure that it truly is a flower and not just a mislabeled grass.


    When there are only two classes (binary classification), the confusion matrix of course would have only four entries. There are special terms that apply to this case when the two classes are “Positive” and “Negative.” Think of diagnosing a disease. Either the patient has the disease or they do not. A doctor can perform a test, and the hope is that the test will determine for sure whether the patient has that disease. Unfortunately, no test is 100% accurate. The two cases in which the test fails are called false positive and false negative. A false positive (also called type I error) occurs when the true state is negative, but the model or test predicts positive. A false negative (also called type II error) is just the opposite: the state is positive, but the model or test predicts a negative. Both errors can be dangerous in fields such as medicine. Some terminology for the four possibilities are shown.


    #figure(table(
  columns: 3, align: left, inset: 6pt,
  table.header([], [#strong[Predicted Positive]], [#strong[Predicted Negative]]),
  [Actual Positive], [True Positive (TP) – hit], [False Negative (FN) – miss],
  [Actual Negative], [False Positive (FP) – false alarm], [True Negative (TN) – correctly rejected],
))


    Adjustments to a model that reduce one type of error generally increases the rate of the other type of error. That is, if you want your model to have a very low rate of false negatives, then it may become positively biased (more sensitive), predicting more positives whether they are true or false positives. On the other hand, if the rate of false positives needs to be reduced, then the model may become negative biased (less sensitive) and yield more negative predictions overall, both true and false negatives.


    
      
==== Visualizing Confusion Matrices


      Confusion matrices are often shaded or colored in a way to show contrast of high and low values, which is helpful in locating any abnormal behavior in the model, often called a #strong[heatmap].  provides a heatmap for the flower/tree/grass example.


      #figure(figph[A heatmap with 9 cells. From left to right, columns are labeled “flower,” “tree,” “grass.” From top to bottom, rows are labeled “flower,” “tree,” “grass.” The top row reads 23, 3, 9. The middle row reads 2, 32, 0. The bottom row reads 12, 1, 18. Cells are varying shades of gray with 0 being the lightest and 32 being the darkest.], alt: "A heatmap with 9 cells. From left to right, columns are labeled “flower,” “tree,” “grass.” From top to bottom, rows are labeled “flower,” “tree,” “grass.” The top row reads 23, 3, 9. The middle row reads 2, 32, 0. The bottom row reads 12, 1, 18. Cells are varying shades of gray with 0 being the lightest and 32 being the darkest.", caption: [Heatmap for Flower/Tree/Grass Example])


      The darker shades indicate higher values. The main diagonal stands out with darker cells. This is to be expected if our classifier is doing its job. However, darker shades in cells that are off the main diagonal indicate higher rates of misclassification. The numbers 12 in the lower left and 9 in the upper right are prominent. The model clearly has more trouble with flowers and grasses than it does with trees.


    
    
      
==== Generating a Confusion Matrix in Python


      The dataset #link("https://openstax.org/r/CollegeCompletionData")[#strong[CollegeCompletionData.csv]] contains 62 data points (from which only 12 points of data were used in ). The following code produces a logistic regression based on all 62 points. Then, the confusion matrix is generated. Finally, a visual display of the confusion matrix as a heatmap is generated using another visualization library called #link("https://openstax.org/r/seaborn")[#strong[seaborn]]. If you do not need the heatmap, then just type print(cf) directly after finding the confusion matrix to print a text version of the matrix.


      #notebox("Note", rgb("#8a94a6"), rgb("#556666"), rgb("#f7f8fa"))[
\# Import libraries
      import pandas as pd  \#\# for dataset management
      import matplotlib.pyplot as plt \#\# for data visualization
      import seaborn as sns \#\# for heatmap visualization
      from sklearn.linear\_model import LogisticRegression
      from sklearn.metrics import confusion\_matrix
      
      \# Read data
      data = pd.read\_csv('CollegeCompletionData.csv').dropna()
      x = data\[\['GPA'\]\]
      y = data\['Completion'\]
      
      \# Build the logistic model
      model = LogisticRegression()
      model.fit(x,y)
      
      \# Generate model predictions
      y\_pred = model.predict(x)
      
      \# Generate the confusion matrix
      cf = confusion\_matrix(y, y\_pred)
      
      \# Plot the heatmap using seaborn and matplotlib
      sns.heatmap(cf, annot=True, fmt='d', cmap='Blues', cbar=True)
      plt.xlabel('Predicted')
      plt.ylabel('True')
      plt.yticks(rotation=0)
      plt.title('Confusion Matrix')
      plt.show()
    
        The resulting output will look like this:


        #figure(figph[A confusion matrix visualizing the performance of a binary classification model. The matrix has two rows and two columns, representing the true and predicted classes, respectively. The diagonal cells (top-left and bottom-right) show the number of correctly classified instances (true positives and true negatives), while the off-diagonal cells show the number of misclassified instances (false positives and false negatives). Clockwise from top left the boxes read: 20, 10, 24, 8. A color scale on the right indicates the frequency of each value. It runs up the right side of the matrix starting as light blue at 5.0 and getting darker at the top at 25. The color of the boxes aligns with this scale.], alt: "A confusion matrix visualizing the performance of a binary classification model. The matrix has two rows and two columns, representing the true and predicted classes, respectively. The diagonal cells (top-left and bottom-right) show the number of correctly classified instances (true positives and true negatives), while the off-diagonal cells show the number of misclassified instances (false positives and false negatives). Clockwise from top left the boxes read: 20, 10, 24, 8. A color scale on the right indicates the frequency of each value. It runs up the right side of the matrix starting as light blue at 5.0 and getting darker at the top at 25. The color of the boxes aligns with this scale.", caption: none)
]


      From the preceding analysis, we find 20 true positives and 24 true negatives, representing an accuracy of #math.equation(block: false, alt: "44 / 62 equals 71 %")[$44 / 62 = 71 %$]. There were 10 false negatives and 8 false positives.