Linear regression
clear, clc
data = [ 10 20 30 40 50 60 70 80
25 70 380 550 610 1220 830 1450];
x = data(1,:);
y = data(2,:);
plot(x,y,'or')
n = length(x); %number of data point
sum(x.*y) %each of "x" multiply by each of "y"
a1 = ( n*sum(x.*y) - sum(x)*sum(y))/( n*sum(x.*x)-sum(x)^2)
ybar = sum(y)/n; xbar =sum(x)/n;
a0 = ybar - a1*xbar
yr = a0 +a1*x;
hold on
plot(x,yr,'b-')
st = sum((y-ybar).^2); % square error wrt. ybar
sr = sum((y-yr).^2); % aquare error wrt. ymodel
rr = (st-sr)/st % r square
เมื่อ RUN แล้วในหน้า command window จะแสดง
ans =
312850
a1 =
19.4702
a0 =
-234.2857
rr =
0.8805
ทำให้เป็น function
function [a0,a1,rr] = fregress_lin(x,y)
n = length(x); %number of data point
a1 = ( n*sum(x.*y) - sum(x)*sum(y))/( n*sum(x.*x)-sum(x)^2);
ybar = sum(y)/n; xbar =sum(x)/n;
a0 = ybar - a1*xbar;
yr = a0 +a1*x;
st = sum((y-ybar).^2); % square error wrt. ybar
sr = sum((y-yr).^2); % aquare error wrt. ymodel
rr = (st-sr)/st ; % r square
ใช้ function ที่ทำไว้มาแก้โจทย์
clear, clc
data = [ 10 20 30 40 50 60 70 80
25 70 380 550 610 1220 830 1450];
x = data(1,:);
y = data(2,:);
plot(x,y,'or')
[a0,a1,rr] = fregress_lin(x,y);
yr = a0 +a1*x;
hold on
plot(x,yr,'b-')
y10 = a0 + a1*10
เมื่อ RUN แล้วในหน้า command window จะแสดง
y10 =
-39.5833
และกราฟ
โจทย์ Enzyme Kinetics
clear,clc
v = [ 1.3 1.8 3 4.5 6 8 9];
s = [ 0.07 0.136 0.22 0.275 0.335 0.35 0.36];
[a0,a1,rr] = fregress_lin(s,v)
vr = a0 + a1*s;
plot (s,v,'or',s,vr,'b-')
เมื่อ RUN แล้วในหน้า command window จะแสดง
a0 =
-1.4240
a1 =
24.9529
rr =
0.8808
และแสดงกราฟ
clear,clc
v = [ 1.3 1.8 3 4.5 6 8 9];
s = [ 0.07 0.136 0.22 0.275 0.335 0.35 0.36];
[a0,a1,rr] = fregress_lin(s,v)
vr = a0 + a1*s;
plot (s,v,'or',s,vr,'b-')
%Michaelis Menten model
y = 1./v;
x = 1./s;
[a0M, a1M, rr] = fregress_lin(x,y)
vm = 1/a0M; ks = vm*a1M;
v_Mic = vm.*s./(ks+s)
plot(s,v,'or',s,vr,'b-',s,v_Mic,'g--')
เมื่อ RUN แล้วในหน้า command window จะแสดง
a0 =
-1.4240
a1 =
24.9529
rr =
0.8808
a0M =
0.0133
a1M =
0.0570
rr =
0.9224
v_Mic =
1.2090 2.3138 3.6729 4.5355 5.4532 5.6789 5.8285
และกราฟ
ในโปรแกรม Matlab มี ฟังก์ชั่นการใช้ curve fitting อยู่แล้ว โดยการ พิมพ์ cftool ในหน้า command window จะแสดงหน้าต่าง
Interpolation
Polynomial interpolation
Example : Determining polynomial coefficient
clear,clc
%data
T = [ 300 400 500];
p = [ 0.616 0.525 0.457];
A = [T(1) 1
T(2) 2];
b = [p(1)
p(2)];
a = inv(A)*b;
TT = 350; %if T = 350 degree
pp = a(1)*TT + a(2)
เมื่อ RUN แล้วในหน้า command window จะแสดง
pp =
0.5705
เปรียบเทียบ 1st กับ 2nd order
clear,clc
%data
T = [ 300 400 500];
p = [ 0.616 0.525 0.457];
% 1st order
A = [T(1) 1
T(2) 1];
b = [p(1)
p(2)];
a = inv(A)*b;
TT = 350;
pp = a(1)*TT + a(2)
%compare to 2nd order
A = [T(1)^2 T(1) 1
T(2)^2 T(2) 1
T(3)^2 T(3) 1];
b = [p(1)
p(2)
p(3)];
a = inv(A)*b;
TT=350;
pp = a(1)*TT^2 + a(2)*TT + a(3)
เมื่อ RUN แล้วในหน้า command window จะแสดง
pp =
0.5705 (1 order)
pp =
0.5676 (2 order)
นำไปทำเป็น function
function [a1,a0,yi] = finterp_lin(x,y,xi)
A = [x(1) 1
x(2) 1];
b = [y(1)
y(2)];
a = inv(A)*b;
a1 = a(1); a0=a(2);
yi = a(1)*xi + a(2);
ใช้ function ที่ทำไว้มาแก้โจทย์
clear,clc
%data
T = [ 300 400 500];
p = [ 0.616 0.525 0.457];
%find P at T=350
[a0,a1,p] = finterp_lin(T,p,350)
เมื่อ RUN แล้วในหน้า command window จะแสดง
a0 =
-9.1000e-004
a1 =
0.8890
p =
0.5705
clear,clc
%data
T = [ 300 400];
p = [ 0.616 0.525];
%find P at T=350
[a0,a1,pr] = finterp_lin(T,p,350)
T = [300 400 500];
p = [0.616 0.525 0.457];
pr = interp1(T,p,350,'spline')
เมื่อ RUN แล้วในหน้า command window จะแสดง
a0 =
-9.1000e-004
a1 =
0.8890
pr =
0.5705
pr =
0.5676
เสร็จจ้าาาาา!!!!
ไม่มีความคิดเห็น:
แสดงความคิดเห็น