Skip to content
Permalink
d2e428c722
Switch branches/tags

Name already in use

A tag already exists with the provided branch name. Many Git commands accept both tag and branch names, so creating this branch may cause unexpected behavior. Are you sure you want to create this branch?
Go to file
 
 
Cannot retrieve contributors at this time
269 lines (237 sloc) 7.24 KB
# 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
```