#1) Create a least-squares function called least_squares.m that accepts a Z-matrix and dependent variable y as input and returns the vector of best-fit constants, a, the best-fit function evaluated at each point, and the coefficient of determination, r2.
function [a,fx,r2]=least_squares(Z,y)
a = (Z'*Z)\(Z'*y);
fx=Z*a;
Sr = sum((y-Z*a).^2);
r2 = 1-Sr/sum((y-mean(y)).^2);
end
Test your function on the sets of data in script problem_1_data.m and show that the following functions are the best fit lines:
a. y=0.3745+0.98644x+0.84564/x
b. y=-11.4887+7.143817x-1.04121 x^2+0.046676 x^3
c. y=4.0046e^(-1.5x)+2.9213e^(-0.3x)+1.5647e^(-0.05x)
The function returns the coefficients of each of the best fit constants in vector a. #2) Use the Temperature and failure data from the Challenger O-rings in lecture_18 (challenger_oring.csv). Your independent variable is temerature and your dependent variable is failure (1=fail, 0=pass). Create a function called cost_logistic.m that takes the vector a, and independent variable x and dependent variable y. a)
function [J, grad]=cost_logistic(a,x,y)
sigma=@(t) 1./(1+exp(-t));
a0=a(1);
a1=a(2);
t=a0+a1*x;
J = sum(-y.*log(sigma(t))-(1-y).*log(1-sigma(t)));
grad(1)= sum((y-sigma(t)));
grad(2)= sum((y-sigma(t)).*x);
end
b)use the following code to solve for a0 and a1
% Set options for fminunc
options = optimset('GradObj', 'on', 'MaxIter', 400);
% Run fminunc to obtain the optimal theta
% This function will return theta and the cost
[theta, cost] = fminunc(@(a)(cost_logisitc(a, x, y)), [0 0], options);
c. plot the data and the best-fit logistic regression model Data Plot
Best-fit logistic regression model
#3) a.Finish the function boussinesq_lookup.m so that when you enter a force, q, dimensions of rectangular area a, b, and depth, z, it uses a third-order polynomial interpolation of the four closest values of m to determine the stress in the vertical direction
function [sigma_z]=boussinesq_lookup(q,a,b,z)
% function that determines stress under corner of an a by b rectangular platform
% z-meters below the platform. The calculated solutions are in the fmn data
% m=fmn(:,1)
% in column 2, fmn(:,2), n=1.2
% in column 3, fmn(:,3), n=1.4
% in column 4, fmn(:,4), n=1.6
fmn= [0.1,0.02926,0.03007,0.03058
0.2,0.05733,0.05894,0.05994
0.3,0.08323,0.08561,0.08709
0.4,0.10631,0.10941,0.11135
0.5,0.12626,0.13003,0.13241
0.6,0.14309,0.14749,0.15027
0.7,0.15703,0.16199,0.16515
0.8,0.16843,0.17389,0.17739];
m=a/z;
n=b/z;
if (n>1.5)
n=1.6;
elseif (1.3<=n)&&(n<=1.5)
n=1.4;
elseif 1.3<=n
n=1.2;
end
f=ones(1,4);
x=ones(1,4);
%for loop to find the values in the fmn matrix that should be used.
for i=1:4
[~,p]=min(abs(m-fmn(:,1)));
M =fmn(p,1);
fmn(p,1)=0;
if n==1.2
t=fmn(p,2);
elseif n==1.4
t=fmn(p,3);
elseif n==1.6
t=fmn(p,4);
end
x(i)= M;
f(i)=t;
end
%coefficients for the polynomial are found with the b values.
b1=f(1);
b2=(f(2)-f(1))/(x(2)-x(1));
b3=(((f(3)-f(2))/(x(3)-x(2)))-((f(2)-f(1))/(x(2)-x(1)))/(x(3)-x(1)));
b4=(((f(4)-f(3))/(x(4)-x(3)))-((f(3)-f(2))/(x(3)-x(2)))-((f(2)-f(1))/(x(2)-x(1))))/(x(4)-x(1));
%the coefficients are plugged into the third order polynoomial
%interpolation.
f3=b1+(b2*(m-x(1)))+(b3*(m-x(1))*(m-x(2)))+(b4*(m-x(1))*(m-x(2))*(m-x(3)));
%Stress in the vertical direction is equal to force*the third order
%polynomial.
sigma_z=q*f3;
end
This function takes a, b, and z and calculates m and n then finds the closest four points and does a fourth order interpolatation to find the function and then multiplies that by the force to find the stress in the vertical direction.
b.Copy the boussinesq_lookup.m to a file called boussinesq_spline.m and use a cubic spline to interpolate in two dimensions, both m and n, that returns sigma_z
fmn= [0.1,0.02926,0.03007,0.03058
0.2,0.05733,0.05894,0.05994
0.3,0.08323,0.08561,0.08709
0.4,0.10631,0.10941,0.11135
0.5,0.12626,0.13003,0.13241
0.6,0.14309,0.14749,0.15027
0.7,0.15703,0.16199,0.16515
0.8,0.16843,0.17389,0.17739];
fmn1=fmn(:,2:4);
x=[1 2 3];
y=[1 2 3 4 5 6 7 8];
sigma_z=interp2(x, y, fmn1, 0.46, 1.4, 'spline')