Statistical Computing and Bayesian Stats Assignment

This assignment is based on statistical computing and Bayesian stats, covering topics such as multinomial regression, correlation function, MLE computation, and ridge regression. Students are required to load sample data, fit a model, compute coefficients, and plot estimates of posterior distributions. The assignment also includes a case study on train cancellation information, where students need to load the data, compute predictor variables, and visualize the results using plots and histograms.

UNIVERSITY AFFILIATION
FACULTY OR DEPARTMENT
STAT878: STATISTICAL COMPUTING
TITLE:
ASSIGNMENT I
STUDENT NAME:
STUDENT REGISTRATION NUMBER:
PROFESSOR (TUTOR)
DATE OF SUBMISSION

1
QUESTION 1: MLE WITH NEWTON-RAPHSON (21 marks)
Part A: Statistical Model
%% STATS878: STATISTICAL COMPUTING
% Student name
% Student Registration number
%% question 1
% load the sample data
% the import only picks the numerical data on temperature, rainfall, and
% ... cancellation respectively.
trainInfo=xlsread('TrainCancellation.xlsx');
disp(trainInfo)
temp=trainInfo(:,1)
rainfall=trainInfo(:,2)
cancellation=trainInfo(:,3)
weather=[temp rainfall];
days=1:1:365;
figure(1)
plot(days, temp,'r',days,rainfall,'b')
title('Weather patterns in March 2017- Feb 2018')
xlabel('days')
ylabel('Weather Patterns')
legend('Temperature','Rainfall');
%% Defining the model, variables, parameters, probabilistic relationship
% defining a nominal response variable using a categorical array
tc=categorical(cancellation);
%% fitting the multinominal regression model to predict the species
... the weather information; mnrfit treats NaNs at missing values and
... ignores them
[B,dev,stats]=mnrfit(weather,tc,'model','hierarchical');
disp(B) % nominal model for the response category relative risks
stats.p % checking the statistical significance of the model
... coefficient.
stats.se % Requests the standard error of coefficient estimates
%% calculating the 95% confidence limits for the coefficients
LL=stats.beta-1.96*stats.se;
UL=stats.beta+1.96*stats.se;
[LL(:,1) UL(:,1)]

2
figure(2)
plot(days,cancellation)
[R,M,B]=regression(temp,cancellation)
figure (3)
plotregression(temp,cancellation)
[s,c,r,d]=regr(temp,cancellation,1,2)
The solution is,
0 50 100 150 200 250 300 350 400
days
0
10
20
30
40
50
60
70
Weather Patterns
Weather patterns in March 2017- Feb 2018
Temperature
Rainfall

3
0 5 10 15 20 25 30
Target
0
5
10
15
20
25
30
Output ~= 0.0027*Target + -0.038
: R=0.10036
Data
Fit
Y = T
Part B: Correlation function
It is important to note that in regression estimates, in an instance where the errors are normally
distributed, the OLS=MLE which is the numerical solution. The regression model is given as
output ≅ 0.0027∗Target ± 0.038
The correlation can be obtained by
x= A ¿ ( OLS )
Part C: MLE computation
%% PART C: MLE
function b=slogit(cancellation,days,b0)
if nargin==2,maxiter=50;b0=(0,0); end
if nargin==3,b0=(0,0); end
oldb=b0';
dx=[ones(length(x),1),x]; % dx=design matrix by attaching 1 to x

4
for i=1:maxiter
eta=dx*oldb;
t1=exp(eta);
p=t1./(1+t1);
t1=p.*(1-p);
g1=sum(y-p);
g2=sum(x.*(y-p));
h11=sum(t2);
h12 = sum(x .* t2);
h22 = sum((x.^2) .* t2);
g = (g1; g2); %gradient
H = [h11, h12; h12, h22];
inc = H \ g; %or inc = inv(H) * g; inc = H-1g
newb = oldb + inc;
if all(abs(newb – oldb) < 1e-4)
break;
else
oldb = newb;
end
end
b = newb;
QUESTION 2: BAYESIAN STATS: PLOT ESTIMATE OF POSTERIOR (14 marks)
%% QUESTION 2: BAYESIAN POSTERIOR PLOTS
function esti = normb2(x, m)
% Model: x ~ N(mu, sig2), mu ~ N(1, 1), sig2 ~ chisq(1)
% x = observed data; m = size of Monte Carlo sample
x = x(:); % turn x into a column vector, if it isn't already
n = length(x);
% Sample mu & sigma:
mu = 1 + sqrt(1) * randn(m, 1); % i.e., mu + sig*N(0, 1)
sig2 = chi2rnd(1, m, 1); sig = sig2.^(1/2);
w = zeros(m, 1); llik = zeros(m, 1); % pre-allocate memory
for j = 1:m
% Compute jth Monte Carlo log-likelihood:
llik(j) = -n*log(sig(j)) - 1/(2*sig(j)^2)*sum((x-mu(j)).^2);
end
for j = 1:m
w(j) = 1./sum(exp(llik-llik(j))); % jth MC weight
end
% Return weighted sums of sampled (mu,sig2) values:
esti = [sum(w.*mu); sum(w.*sig2)];
k = 2; n = 10; %defines values for k & n
t = 1; rnd = zeros(500, 1); %generate 500 r.n.
while t <= 500
p = rand;
ratio = (p^k * (1-p)^(n-k))/((k/n)^k * (1-k/n)^(n-k));
u = rand;
if u <= ratio
rnd(t) = p;
t = t + 1;

5
end
end
x = 0:0.01:1;
y = x.^(3-1) .* (1-x).^(9-1) ./ beta(3, 9);
histogram(rnd, 0:0.04:1, 'Normalization', 'pdf');
hold on; plot(x, y, 'r'); hold off; % add red curve
The solution is obtained as,
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
0
0.5
1
1.5
2
2.5
3
3.5
4
QUESTION 3: RIDGE REGRESSION (15 marks)
Procedure:
(i) Loading the train cancellation information
(ii) The predictor variables are temperature and rainfall, and the response variable is train
cancellation
(iii) The predictor variables are plotted against each other.
(iv) The correlation between the predictor variables is obtained with respect to each other.
The coefficient estimates for a multilinear model with interaction terms are computed
for a range of ridge parameters.
(v) The plots for ridge regression are obtained as shown in the output plots.
%% question 3

6
% Ridge regression
tc=xlsread('TrainCancellation.xlsx');
temp=tc(:,2);
rainfall=tc(:,3);
cancellation=tc(:,4);
weather=[temp rainfall];
D=x2fx(weather, 'interaction')
D(:,1)=[]; %No constant term
k=1:1:365;
b=ridge(cancellation,D,k)
%% plotting the ridge trace
figure(1)
plot(k,b,'LineWidth',2)
ylim([-100 100])
grid on
xlabel('Ridge Parameter')
ylabel('Standardized Coefficient')
title('{\bf Ridge Trace}')
legend('Temperature','Rainfall','Temperature-Rainfall','Cancellations')
REFERENCES

7
[1] Dobson, A. J., and A. G. Barnett. An Introduction to Generalized Linear Models. Chapman
and Hall/CRC. Taylor & Francis Group, 2008.
[2] K. Sigmon and T. A. Davis, MATLAB Primer, Sixth Edition, Chapman and
Hall/CRC, Boca Raton, FL, 2002
[3] Marquardt, D.W. “Generalized Inverses, Ridge Regression, Biased Linear Estimation, and
Nonlinear Estimation.” Technometrics. Vol. 12, No. 3, 1970, pp. 591–612.