Bayesian regression problem

Best-fit model is not necessarily the best model. It is important to balance between a good fit to the data and model complexity. The purpose of this exercise is to illustrate this idea via a regression problem, which was discussed.

Table 1 (see appendix) contains 8 observations (values stored as x and y vectors in INPUT.mat). Use load ('INPUT.mat') in MATLAB to read in these values. This matlab file can be downloaded from CCLE course website /problem sets.

We can define all possible polynomial regression models as:

y = β₀ + β₁x + β₂x² + ? + β_p x^p + ε, where ε ~ Normal (0, σ²)

In this exercise, we consider seven possible models: p = 0, 1, 2,..., 6.

Our goal is to decide which one of the seven models is the "best" one to explain the observed data. For each model with specific parameters β, "goodness of fit" can be measured by the likelihood term P(y|β,M). Here, β is a vector of regression coefficients, and M indicates a polynomial regression model of order p. Model evidence is evaluated using P(y|M), which describes how likely the data are generated by a polynomial model.

(a) Analyze the goodness of fit of each model M by computing the likelihood for each model, based on the predefined regression coefficients, b provided in INPUT.mat (also listed in Table 2 for each model in appendix, see more detail in appendix).

In statistical terms, likelihood can be understood as how the data are generated/sampled from a model. The assumption that ε follows a normal distribution with zero mean and a constant standard deviation σ tells you the variation in data generation.

We can write the likelihood probability distribution as

y_{i ~} Normal(y ^_i, σ²), where y ^_I = b₀ + b₁x_i + b₂x_i² + ? + b_p x_i^p

Here, y ^_i (called "y-hat") is the predicted value for the i^th observation on y. We further assume that each data point is independently sampled. Then, for a particular regression model M with order p, the likelihood is given by

P(y|β, M)= _i=1∏⁸?(y_i; y ^_i, σ²)

?(x; μ, σ²) refers to the probability density function of the normal distribution (i.e., norm pdf function in MATLAB) with mean μ and standard division σ. Please use σ=5 for your likelihood calculation.

Since likelihoods across different models could differ by orders of magnitude, for a better illustration, it is more advantageous to plot the natural logarithm of the likelihoods instead of the raw likelihood values. Present a plot of log-likelihood against the orders of polynomial p. What trend do you observe from the log-likelihood plot? Which model gives you the "best fit"?

(b) Evaluate each model M by its model evidence P(y|M), which is given by

P(y|M)= -∞∫^+∞P(y|β, M)P(β)dβ

Computing this integral analytically is hard. Instead, we use the discrete approximation:

-∞∫^+∞P(y|β, M)P(β)dβ  ≈  1/N _j=1∑^N P(y|β_j, M)

To simplify your calculation, we assume the prior P(β) to be a uniform distribution, i.e., β_k~Uniform(A,B), where A = b_k - 0.5, and B = b_k + 0.5, for k = 0, 1, ..., p (p is the order of the polynomial regression of model M, b_k is the k^th value in Table 2 for each Model M).

Using sampling approach to implement the Bayesian model. Sample N sets of β values for each model M according to the prior distribution. For each sampled β_j, compute P(y|β_j, M), which is given by the likelihood equation. Use N = 500.

Present a bar chart of model evidence against the orders of polynomial p. Which model gives you the highest model evidence?

(c) Write a short paragraph to discuss which regression model is the best for this set of data.

Attachment:- Assignment.rar