%matlab code attached below:
% if setting x(t) = [1 : m] , then f(t) is the following array
% (the goal is to find the coefficients in the approximation model.)
f=[13, 13, 26, 52, 52, 102, 145, 174, 205, 249, 257, 263, 282, 294, 303, 321, 344, 358, 340, 276, 214, 225, 237, 255, 279, 288, 272, 225, 219, 228, 251, 268, 279, 271, 266, 232, 216, 222, 234, 254, 269, 264, 230, 195, 187, 193, 208, 218, 237, 239, 220, 198, 195, 207, 214, 223, 231, 239, 233, 212, 208, 218, 228, 237, 244, 244, 240, 226, 200, 200, 213, 232, 241, 243, 240, 231, 236, 260, 282, 302, 315, 325, 307, 295, 289, 295, 305, 321, 335, 313, 289, 264, 263, 263, 277, 290, 291, 284, 271, 242, 236, 241, 256, 273, 274, 267, 247, 235, 242, 267, 285, 303, 321, 325, 312, 307, 304, 309, 331, 355, 364, 354, 338, 313, 305, 306, 319, 322, 336, 336, 317, 299, 295, 298, 317, 342, 352, 364, 361, 366, 375, 379, 391, 410, 425, 429, 414, 391, 389, 410, 437, 457, 450, 439, 412, 395, 389, 406, 425, 453, 470, 444, 408, 409, 410, 423, 442, 464, 475, 472, 457, 438, 428, 441, 456, 471, 473, 465, 443, 407, 409, 428, 444, 460, 463, 461, 447, 420, 408, 413, 441, 462, 457, 437, 423, 409, 389, 399, 416, 433, 430, 426, 400, 384, 382, 395, 415, 432, 440, 445, 442, 432, 433, 462, 489, 538, 560, 564, 539, 534, 521, 523, 547, 575, 590, 583, 555, 541, 537, 543, 558, 575, 577, 570, 556, 538, 520, 526, 543, 558, 562, 553, 527, 507, 495, 492, 507, 527, 539, 531, 517, 510, 508, 511, 523, 554, 562, 554, 544, 537, 541, 546, 553, 567, 566, 554, 528, 516, 522, 510, 523, 542, 563, 570, 558, 553, 559, 557, 561, 581, 596, 585, 599, 583, 578, 562, 553, 558, 584, 589, 587, 590, 582, 580, 590, 624, 640, 649, 653, 647, 619, 624, 631, 659, 678, 688, 695, 681, 657, 632, 644, 674, 703, 733, 741, 724, 701, 686, 677, 687, 697, 707, 694, 682, 657, 642, 645, 666, 675, 675, 666, 655, 644, 647, 655, 675, 691, 692, 696, 686, 685, 687, 696, 705, 714, 703, 697, 697, 700, 689, 707, 729, 746, 755, 761, 768, 775, 780, 786, 807, 813, 809, 811, 807, 813, 796, 791, 807, 815, 818, 799, 809, 811, 791, 791, 781, 773, 781, 777, 784, 801, 787, 785, 784, 805, 809, 804, 807, 819, 818, 813, 813, 833, 830, 840, 851, 844, 837, 834, 842, 854, 859, 865, 864, 860, 847, 841, 855, 861, 879, 893, 888, 873, 837, 814, 812, 819, 837, 849, 862, 862, 874, 873, 889, 892, 915, 931, 930, 921, 918, 914, 892, 894, 899, 896, 874, 879, 876, 864, 849, 859, 869, 879, 872, 875, 873, 831, 816, 834, 858, 868, 874, 880, 879, 878, 877, 887, 913, 930, 935, 945, 959, 976, 990, 1003, 1023, 1019, 1003, 993, 985, 993, 998, 971, 958, 971, 984, 994, 1013, 1019, 1020, 1018, 1029, 1025, 1028, 1045, 1042, 1043, 1028, 1027, 1031, 1032, 1042, 1037, 1007, 996, 973, 979, 977, 992, 1006, 1007, 996, 996, 995, 992, 1004, 1013, 1015, 1012, 1011, 995, 1000, 1003, 987, 990, 1001, 975, 946, 949, 949, 952, 955, 964, 958, 972, 967, 953, 960, 959, 970, 967, 964, 959, 934, 927, 915, 919, 930, 958, 986, 997, 998, 991, 981, 972, 982, 990, 992, 991, 995, 975, 973, 977, 997, 1026, 1022, 1006, 1006, 1005, 1000, 1002, 988, 1002, 1018, 1019, 1006, 983, 978, 989, 975, 991, 1016, 1010, 1016, 1019, 1033, 1040, 1036, 1061, 1069, 1067, 1059, 1055, 1056, 1056, 1045, 1060, 1073, 1081, 1062, 1055, 1057, 1040, 1035, 1032, 1041, 1044, 1053, 1034, 1015, 976, 963, 984, 1001, 1018, 1028, 1033, 1025, 1024, 1031, 1045, 1066, 1079, 1095, 1081, 1084, 1116, 1127, 1134, 1158, 1159, 1174, 1158, 1158, 1155, 1156, 1163, 1190, 1180, 1164, 1142, 1146, 1126, 1113, 1128, 1151, 1159, 1136, 1131, 1115, 1096, 1079, 1079, 1086, 1094, 1075, 1077, 1063, 1068, 1047, 1054, 1073, 1096, 1094, 1078, 1066, 1057, 1054, 1049, 1066, 1085, 1091, 1078, 1077, 1067, 1087, 1104, 1115, 1142, 1139, 1143, 1155, 1151, 1157, 1160, 1147, 1165, 1164, 1148, 1131, 1113, 1105, 1105, 1098, 1117, 1131, 1132, 1145, 1137, 1127, 1128, 1123, 1130, 1135, 1152];
m = length(f);
f = reshape(f, m, 1);
x = [1 : m]';
%%
%% find linear approximation: f1(t) = c(1) + c(2)*x(t)
%%
%% (the following single line contains the essense of a linear least square method.
%% you can easily change it into a script function)
A = [ones(m,1), x]; c = A \ f;
f1 = c(1) + c(2).* x;
%comute the relative root-mean-square (RMS) error:
RMS_linear_error = sqrt(sum((f-f1).^2)/m)/norm(f);
fprintf(' Linear Model relative RMS error = %e\n', RMS_linear_error);
%% for this specific x(t)=t, the f(t) is simpplified into
%% f(t) = - a/N * tˆ2 + a* t - 1
%% but this model does not lead to a useful least-square solution.
%% this is because the linear approximation: f2(t) = c(1) t + c(2)* tˆ2 - 1
%% would lead to a matrix of just rank-1, thus not good for a least square solution.
%% modify model so that the three terms can be used to generate three coefficients (instead of just two),
%% this turns the approximation function into a quadractic model:
%%
%% the goal is to find the quadratic approximation: f3(t) = c(1) + c(2)*x(t) + c(3)*xˆ2
%%
%% equivalently (trying to be close to the original model),
%% use the model f(t)=-a*x(t)^2/N +(a+1)*x(t) - b*x(t+1), trying to find a, N, b.
%%
%% i.e., solve for f(t) = c(3) t^2 + c(2) t + c(1) = -a/N t^2 + (a -b +1) *t - b ,
%%
%% this always has a least square solution, then one can solve for the a, N, and b by solving
%% -a/N = c(3), a - b + 1 = c(2), -b = c(1)
%%
%% (the following single line contains the essense of a quadractic least square method,
%% you can easily change it into a script function)
A = [ones(m,1), x, x.*x]; c = A \ f
f3 = c(1) + c(2).* x + c(3).*x.*x ;
RMS_quadratic_error = sqrt(sum((f-f3).^2)/m)/norm(f);
fprintf(' Quadratic Model relative RMS error = %e\n', RMS_quadratic_error);
%
% draw the two approximations
%
plot(x, f, 'b:', x, f1, 'r--', x, f3, 'g-')
legend('data f(t)', ['linear, rRMS=', num2str(RMS_linear_error)], ...
['quadratic, rRMS=', num2str(RMS_quadratic_error)], 'Location', 'Best')
print(gcf, '-dpng', 'ntest.png')
% screen output the a, b, N
b = - c(1)
a = c(2) - 1 + b
N = -a/c(3)
扫码加好友，拉您进群
全部版块
我的主页
收藏
