Probit regression is a machine learning technique for binary classification problems. For example, the Banknote authentication problem is to predict if a banknote (think dollar bill or euro) is real (0) or a forgery (1), based on four predictor values of a digital image of the banknote (variance, skewness, kurtosis, entropy).

Probit regression is very similar to logistic regression. In logistic regression, you feed numeric predictor values as input. The output value is between 0.0 and 1.0 where an output of less than 0.5 indicates class 0 and an output of greater than or equal to 0.5 indicates class 1. In logistic regression, the output is computed using a very simple function called logistic_sigmoid().

The phi() function is very similar to the logistic_sigmoid() function.

Probit regression works exactly like logistic regression except that instead of using the logistic_sigmoid() function, you use the phi() function. The phi() function computes the area under the curve of the standard Normal (aka Gaussian) distribution.

I coded up a quick demo. I used the Banknote Dataset where I normalized the four predictor values by dividing each by 20. I split the data into a 1,000-item training set and a 372-item test set. The raw Banknote data can be found at https://archive.ics.uci.edu/ml/datasets/banknote+authentication. The normalized data looks like:

# variance  skewness   kurtosis   entropy  class
-0.177550   0.094775   0.009325  -0.122045   1
 0.065570   0.227310   0.114675   0.011271   0
-0.200865  -0.415615   0.622735  -0.071875   1
-0.255950   0.332430  -0.002499  -0.326030   1
 0.181445   0.040661   0.081385   0.038814   0
. . .

The phi() function is tricky to compute. I used an equation called A and S 7.1.26 (“Abramowitz and Stegun”). The input is any value, the output is between 0.0 and 1.0.

To train the model I used a variation of basic stochastic gradient descent. SGD depends on the gradient of the phi() function which is quite complicated so instead I used SGD for logistic_sigmoid(). My idea is that the two functions are so similar, the gradient for phi() will be very close to the gradient for sigmoid().

Good fun.

The term “probit” in probit regression means “probability unit”. Here are three computer generated images that rely on probability.

Demo code. Replace “lt”, “gt”, “lte”, “gte” with Boolean operator symbols.

# banknote_probit_reg.py

# predict real (0) or forgery (1) from
# variance, skewness, kurtosis, entropy (all div by 20.0)
# data:
# -0.177550  0.094775  0.009325  -0.122045   1
#  0.065570  0.227310  0.114675   0.011271   0

# Anaconda3-2020.02  Python 3.7.6
# Windows 10/11

import numpy as np

# -----------------------------------------------------------

def compute_output(w, b, x):
  # input x, using weights w and bias b
  z = 0.0
  for i in range(len(w)):
    z += w[i] * x[i]
  z += b
  # p = logistic_sigmoid(z)  # logistic regression
  p = phi(z)                 # probit regression
  return p

# -----------------------------------------------------------

def phi(z):
  # cumulative density function for standard Gaussian
  # A_and_S 7.1.26 erf()

  if z "lt" -4.0:
    return 0.0
  if z "gt" 4.0:
    return 1.0

  a0 = 0.3275911
  a1 = 0.254829592
  a2 = -0.284496736
  a3 = 1.421413741
  a4 = -1.453152027
  a5 = 1.061405429

  sign = 0
  if z "lt" 0.0:
    sign = -1;
  else:
    sign = 1;

  x = np.abs(z) / np.sqrt(2.0)
  t = 1.0 / (1.0 + a0 * x);
  erf = 1.0 - (((((a5 * t + a4) * t) + a3) * t + a2) * t + a1) \
    * t * np.exp(-x * x)
  return 0.5 * (1.0 + (sign * erf))

# -----------------------------------------------------------

def logistic_sigmoid(z):
  if z "lt" -4.0:
    return 0.0
  elif z "gt" 4.0:
    return 1.0
  else:
    return 1.0 / (1.0 + np.exp(-z))

# -----------------------------------------------------------

def accuracy(w, b, data_x, data_y):
  n_correct = 0; n_wrong = 0
  for i in range(len(data_x)):
    x = data_x[i]  # inputs
    y = int(data_y[i])  # target 0 or 1
    p = compute_output(w, b, x)
    if (y == 0 and p "lt" 0.5) or (y == 1 and p "gte" 0.5):
      n_correct += 1
    else:
      n_wrong += 1
  acc = (n_correct * 1.0) / (n_correct + n_wrong)
  return acc

# -----------------------------------------------------------

def mse_loss(w, b, data_x, data_y):
  sum = 0.0
  for i in range(len(data_x)):
    x = data_x[i]  # inputs
    y = int(data_y[i])  # target 0 or 1
    p = compute_output(w, b, x)
    sum += (y - p) * (y - p)
  mse = sum / len(data_x)
  return mse

# -----------------------------------------------------------

def main():
  # 0. get ready
  print("\nBegin probit regression with raw Python demo ")
  np.random.seed(1)

  # 1. load data
  print("\nLoading Banknote train and test data to memory ")
  # variance, skewness, kurtosis, entropy (all div by 20.0)
  # 0 = real, 1 = forgery
  # -0.177550  0.094775  0.009325  -0.122045   1
  #  0.065570  0.227310  0.114675   0.011271   0

  train_file = ".\\Data\\banknote_train.txt"
  train_xy = np.loadtxt(train_file, usecols=range(0,5),
    delimiter="\t", comments="#",  dtype=np.float32) 
  train_x = train_xy[:,0:4]
  train_y = train_xy[:,4]

  test_file = ".\\Data\\banknote_test.txt"
  test_xy = np.loadtxt(test_file, usecols=range(0,5),
    delimiter="\t", comments="#", dtype=np.float32)
  test_x = test_xy[:,0:4]
  test_y = test_xy[:,4]

# -----------------------------------------------------------

  # 2. create model
  print("\nCreating probit regression model ")
  wts = np.zeros(4)  # one wt per predictor
  lo = -0.01; hi = 0.01
  for i in range(len(wts)):
    wts[i] = (hi - lo) * np.random.random() + lo
  bias = 0.00

# -----------------------------------------------------------

  # 3. train model
  lrn_rate = 0.01
  max_epochs = 100
  indices = np.arange(len(train_x))  # [0, 1, .. 999]
  print("\nTraining using SGD with lrn_rate = %0.4f " % lrn_rate)
  for epoch in range(max_epochs):
    np.random.shuffle(indices)
    for i in indices:
      x = train_x[i]  # inputs
      y = train_y[i]  # target 0.0 or 1.0
      p = compute_output(wts, bias, x)

      # update all wts and the bias
      for j in range(len(wts)):
        wts[j] += lrn_rate * x[j] * (y - p)  # target - oupt
      bias += lrn_rate * (y - p)
    if epoch % 10 == 0:
      loss = mse_loss(wts, bias, train_x, train_y)
      print("epoch = %5d  |  loss = %9.4f " % (epoch, loss))
  print("Done ")

# -----------------------------------------------------------

  # 4. evaluate model
  print("\nEvaluating trained model ")
  acc_train = accuracy(wts, bias, train_x, train_y)
  print("Accuracy on train data: %0.4f " % acc_train)
  acc_test = accuracy(wts, bias, test_x, test_y)
  print("Accuracy on test data: %0.4f " % acc_test)

  # 5. use model
  print("\nPrediction for dummy [0.2, 0.3, 0.5, 0.7] banknote: ")
  x = np.array([0.2, 0.3, 0.5, 0.7], dtype=np.float32)
  p = compute_output(wts, bias, x)
  print("%0.8f " % p)
  if p "lt" 0.5:
    print("class 0 (real) ")
  else:
    print("class 1 (forgery) ") 

  # 6. TODO: save trained weights and bias to file

  print("\nEnd Banknote probit regression demo ")

if __name__ == "__main__":
  main()