Question

What is the real-life example of the benefit and application of the benefit of Bayesian regression?

Having read the items and it looks having the range of inference (possible values and likelihood) available is the benefit.

• 
  What's the advantages of bayesian version of linear regression, logistic regression etc (1)


• 
  compare bayesian linear regression VS linear regression [closed] (2)


• 
  第14回　ベイズ線形回帰を実装してみよう (implement bayesian regression) (3)

But how it will be utilised in real life (finance, engineering, ...)? I suppose without real-life applications of the benefit, it would not be so useful.

It seems both linear regression and Bayesian regression can produce similar predictions as below.

According to 3, the predictive distribution can give the confidence on the prediction if it is within the dense-color area because of the data is dense, but not in sparse area, eg. prediction at x=5 may not be trustworthy. But I suppose there should be more than that.

What does "having statistical inference" make a difference and how it is utilised in real life? Or Is it a matter of choice to use linear or Bayesian?

Code

Taken from (3) and made a few changes. I am not completely familiar with the logic so if there are mistakes, please correct.

import math

import numpy as np

import matplotlib.pyplot as plt



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

# Gaussian basis function

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

def gaussian(mean, sigma):

    """

    Args:

        mean:  

        sigma: 

    """

    def _gaussian(x):

        return np.exp(-0.5 * ((x - mean) / sigma) ** 2)

    return _gaussian



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

# Design matrix

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

def phi(f, x): 

    bias = np.array([1])  # bias parameter(i = 0)

    # bias+basis

    return np.append(bias, f(x))



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

# Data generation utility

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

from numpy.random import rand

def uniform_variable_generator(samples):

    _random = rand(samples)

    return _random



def noise_generator(samples, mu=0.0, beta=0.1):

    noise = np.random.normal(mu, beta, samples)

    return noise



def sigmoid(x):

  return 1 / (1 + np.exp(-x))



def generator_t_function(x):

    #return np.sin(x)

    return sigmoid(x)



def generator_X_function(x):

    return 2 * np.pi * x

    #return 2 * np.pi * x



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

# Observations

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

#X = np.array([0.02, 0.12, 0.19, 0.27, 0.42, 0.51, 0.64, 0.84, 0.88, 0.99])

#t = np.array([0.05, 0.87, 0.94, 0.92, 0.54, -0.11, -0.78, -0.79, -0.89, -0.04])

samples = 20



#X = np.array([0.02, 0.12, 0.19, 0.27, 0.42, 0.51, 0.64, 0.84, 0.88, 0.99])

#t = np.array([0.05, 0.87, 0.94, 0.92, 0.54, -0.11, -0.78, -0.79, -0.89, -0.04])

X = generator_X_function(uniform_variable_generator(samples))

t = generator_t_function(X) + noise_generator(samples, beta=0.1)



MAX_X = max(X)

NUM_X = len(X)

MAX_T = max(t)

NUM_T = len(t)



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

# Gaussian basis function parameters

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

sigma = 0.1 * MAX_X



# mean of gaussian basis function (11 dimension w1, w2, ... w11)

mean = np.arange(0, MAX_X + sigma, sigma)  



# Basis function

f = gaussian(mean, sigma)



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

# Design matrix

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

PHI = np.array([phi(f, x) for x in X])



#alpha = 0.1 

#beta = 9.0 

alpha = 0.5  # larger alpha gives smaller w preventing overfitting (0 -> same with linear regression)

beta = 5   # Small beta allows more variance (deviation)



Sigma_N = np.linalg.inv(alpha * np.identity(PHI.shape[1]) + beta * np.dot(PHI.T, PHI))



mean_N = beta * np.dot(Sigma_N, np.dot(PHI.T, t))



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

# Bayesian regression

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

xlist = np.arange(0, MAX_X, 0.01) 

plt.title("Bayesian regression")

plt.plot(xlist, [np.dot(mean_N, phi(f, x)) for x in xlist], 'b') 

plt.plot(X, t, 'o', color='r') 

plt.show()



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

# Linear regression

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

# w for linear regression parameter

#w = np.linalg.solve(np.dot(PHI.T, PHI), np.dot(PHI.T, t))

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

l = 0.05

regularization = np.identity(PHI.shape[1])

w = np.linalg.solve(

    np.dot(PHI.T, PHI) + (l * regularization), 

    np.dot(PHI.T, t)

)



xlist = np.arange(0, MAX_X, 0.01) 

plt.title("Linear regression")

plt.plot(xlist, [np.dot(w, phi(f, x)) for x in xlist], 'g') 

plt.plot(X, t, 'o', color='r') 

plt.show()



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

# Predictive Distribution

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

def normal_dist_pdf(x, mean, var): 

    return np.exp(-(x-mean) ** 2 / (2 * var)) / np.sqrt(2 * np.pi * var)



def quad_form(A, x):

    return np.dot(x, np.dot(A, x))



xlist = np.arange(0, MAX_X, 0.01)

#tlist = np.arange(-1.5 * MAX_T, 1.5 * MAX_T, 0.01)

tlist = np.arange(

    np.mean(t) - (np.max(t)-np.min(t)), 

    np.mean(t) + (np.max(t)-np.min(t)), 

    0.01

)

z = np.array([

    normal_dist_pdf(tlist, np.dot(mean_N, phi(f, x)),

    1 / beta + quad_form(Sigma_N, phi(f, x))) for x in xlist

]).T



plt.contourf(xlist, tlist, z, 5, cmap=plt.cm.coolwarm)

plt.title("Predictive distribution")

plt.plot(xlist, [np.dot(mean_N, phi(f, x)) for x in xlist], 'r')

plt.plot(X, t, 'go')

plt.show()

We use a Bayesian model in forecasting retail sales with the SAP Unified Demand Forecast (UDF). The Bayesian approach offers two advantages:

1. We can use priors. For instance, suppose your supermarkets stocks a new item, and we would like to forecast its first Christmas sales. We can simply use the average effect of similar items as a prior for the new item's Christmas effect. Which has the added benefit that the prior gets automatically updated once we have seen the new item's first Christmas sales, so both the prior and the item's own data influence the prediction for the second Christmas.




2. The priors regularize. Our model is heavily over-parameterized, with seasonality, day of week, trend, holidays and tons of promotion predictors, so regularization is hugely important to keep the predictions under control. (Note that the pictures you show in your question do not show such an overparameterized model.)

Aspect 1 is more important to our marketers and to business users, and aspect 2 is more important to me and to statistician or data scientist users.

Yes, there are ways to address both issues without a Bayesian model, e.g., by pooling for aspect 1, and the Lasso or Elastic Net for aspect 2. We just chose the Bayesian approach, and I find the updating particularly elegant.

So, yes, Bayesian regression is used in real life, and it has paid my salary for the last couple of years, along with the salaries of a number of my colleagues. And while making sure your supermarket does not run out of your favorite shampoo.