curve_fitting
1
linear_squares.m
function [ a,fx,r_2 ] = linear_squares(Z,y)
%linear-squares accepts a Z-matrix and dependent variable y and returns the
%vector of best-fit constants. This function allows for more reliable modeling of
%non-linear functions.
a = (Z'*Z)\(Z'*y); %calculate the parameters of function
fx = Z*a;
Sr = sum((y-Z*a).^2); %sum of the squares of the residuals
r_2 = 1-Sr/sum((y-mean(y)).^2); %calculate each elements distance from the regressionn line
end
1a.
%----Calculation of 1a.----------------------------
x = [1 2 3 4 5]'; %Column vector x
y = [2.2 2.8 3.6 4.5 5.5]'; % Column vector y
Z = [x.^0 x x.^-1]; %Compute Z matrix (from book)
[a, fx, r_2] = linear_squares(Z,y) %apply to linear_squares.m
x_fn = linspace(min(x),max(x));
% Plot to show fit
plot(x,y,'o',x_fn, a(1)+a(2)*x_fn+a(3)*x_fn.^-1);
title('Best fit case A');
xlabel('x');
ylabel('y');
Results:
a =
0.3745
0.9864
0.8456
fx =
2.2066
2.7702
3.6157
4.5317
5.4758
r_2 =
0.9996
![fit1a] (fit1a.png)
1b.
%----Calculation of 1b.----------------------------
%Z matrix calculation altered to best fit polynomial
x = [3 4 5 7 8 9 11 12]';
y = [1.6 3.6 4.4 3.4 2.2 2.8 3.8 4.6]';
Z = [x.^0 x.^1 x.^2 x.^3]; %Compute Z matrix (from book)
[a, fx, r_2] = linear_squares(Z,y) %apply to linear_squares.m
x_fn = linspace(min(x),max(x));
plot(x,y,'o', x_fn, a(1)+a(2)*x_fn.^1+a(3)*x_fn.^2+a(4)*x_fn.^3);
title('Best fit case B');
xlabel('x');
ylabel('y');
Results:
a =
-11.4887
7.1438
-1.0412
0.0467
fx =
1.8321
3.4145
4.0347
3.5087
2.9227
2.4947
3.2330
4.9595
r_2 =
0.8290
![fit1b] (fit1b.png)
1c.
%----Calculation of 1c.----------------------------
%Z matrix calculation altered to best fit polynomial
x = [0.5 1 2 3 4 5 6 7 9]';
y = [6 4.4 3.2 2.7 2.2 1.9 1.7 1.4 1.1]';
Z = [exp(-1.5.*x) exp(-0.3.*x) exp(-0.05.*x)]; %Compute Z matrix (from book)
[a, fx, r_2] = linear_squares(Z,y) %apply to linear_squares.m
x_fn = linspace(min(x),max(x));
% Plot to show fit
plot(x,y,'o',x_fn, a(1)*exp(-1.5.*x_fn)+a(2)*exp(-0.3.*x_fn)+a(3)*exp(-0.5.*x_fn));
title('Best fit case C');
xlabel('x');
ylabel('y');
Results:
a =
4.0046
2.9213
1.5647
fx =
5.9321
4.5461
3.2184
2.5789
2.1709
1.8726
1.6425
1.4605
1.1940
r_2 =
0.9971
![fit1c] (fit1c.png)
2
2a
Function
function [cost, gradient] = cost_logisticV3(a, x, y)
% cost_logisticV3 Compute cost and gradient for logistic regression
% cost = cost_logistic(theta, X, y) computes the cost of using theta as the
% parameter for logistic regression and the gradient of the cost
% w.r.t. to the parameters.
% ====================== YOUR CODE HERE ======================
% Compute the cost of a particular choice of a. Compute the partial derivatives and set grad to the partial derivatives of the cost w.r.t. each parameter in theta
cost = 0; %initialize values
gradient = 0; %initalize a gradient
t = a(1).*(x.^0)+a(2).*x; %time function
sigma = 1./(1+exp(-t)); %calculate sigma
cost = sum(-y.*log(sigma)- (1-y).*log(1-sigma)); %Calculation of cost
costFunction = @ (a) sum(-y.*log((1./(1+exp(-(a(1)+a(2).*x)))))-(1-y).*log(1-(1./(1+exp(-(a(1)+a(2).*x))))));
gradient = (1/length(x))*sum((sigma-y).*t); %Calculation of gradient
ai = [0 0]; %initialize a value for a
%---------------------- Plot of Regression -----------------------------
plot(x,y,'x', x, sigma);
title('Regression')
xlabel('Temp (Degrees Fahrenheit)')
ylabel('Pass or Fail (1 or 0)')
end
Script to optimize coefficients of A
%--------- Challenger_cost_logistic.m ---------------------------
%The purpose of the script is to accept the O-ring data, and run it through
%a cost analysis that will also output a plot of the best fit model for
%part failures.
%--- Compile data from Challenger O-ring file -----------
challenger = [53,1;57,1;58,1;63,1;66,0;66.8,0;67,0;67.2,0;68,0;69,0;69.8,1;69.8,0;70.2,0;70.2,0;72,0;73,0;75,0;75,1;75.8,0;76.2,0;78,0;79,0;81,0];
x = challenger(:,1); % Independent variable x assigned to temperature
y = challenger(:,2); % Dependent variable y assigned to status
% 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_logisticV3(a, x,y)), [20 -3]); %Output the cost and gradient of data and graph and values of theta.
2b
theta = [a0,a1] = [16.9868,-0.2659]
2c
![Price of Tea in China] (Price of Tea in China.png)
3
3a
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(:,2), n=1.4
% in column 4, fmn(:,2), n=1.6
fmn= [0.1,0.02926,0.03007,0.03058 %provided data
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; %backslash operator to solve m
n=b/z; %backslash for n
%------------- Building indices ---------------------------------------
%Two loops designed to isolate the idices in the matrix that will be
%interpolated between in the future
mindex = [];
b = sort(abs(fmn(:,1)-m));
close2m = b(1:4);
%Build indices for m values
for j=1:4
for i=1:8
if close2m(j) == abs(fmn(i,1)-m)
mindex=[mindex;i];
end
end
end
%Build indices for desired n values
if n<1.2
n=2;
elseif n>1.6
n=4;
else
n=int8((round(n*5)/5-1)/0.2+1);
end
%----------Value organization --------------------------------------
%reference values of the previously found indices for n & m
c=zeros(4,1);
close2mtrue=zeros(4,1);
for i=1:4 %organize values from the matrix
close2mtrue(i) = fmn(mindex(i),1);
c(i) = fmn(mindex(i),n);
end
sigma_z = q*ppval(interp1(close2mtrue,c,'cubic','pp'),m); %interpolation between the closest n & m values
end
3b
function sigma_z=boussinesq_splineV2(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(:,2), n=1.4
% in column 4, fmn(:,2), n=1.6
fmn= [0.1,0.02926,0.03007,0.03058 %provided data
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; %backslash operator to solve m
n=b/z; %backslash for n
m1 = ppval(spline(fmn(:,1),fmn(:,2)),m); % int. wrt to column 2
m2 = ppval(spline(fmn(:,1),fmn(:,3)),m); % int. wrt to column 3
m3 = ppval(spline(fmn(:,1),fmn(:,4)),m); % interpolate wrt column 4
sigma_z = q*ppval(spline([1.2;1.4;1.6],[m1;m2;m3]),n); %interpolates across all previous interpolations
end