1. % Plot function
2. function [] = draw(x, x_plot, y, f, p_y, x_pred, y_pred, y_true)  
3.     figure;
4.     plot(x_plot, f, 'b');
5.     hold on;
6.     plot(x, y, 'ro');
7.     hold on;
8.     plot(x_pred, y_pred, 'bo', 'LineWidth', 1.5);
9.     hold on;
10.     plot(x_pred, y_true, 'ro', 'LineWidth', 1.5);
11.     hold on;
12.     plot(x, p_y, 'bo');
13.     title('Curve fitting');
14.     xlabel('Days');
15.     ylabel('Drug released (ug)');
16. end
17.  
18. % Residual calculate function
19. function [] = res(x, p_y, y)
20.     residual = zeros(1, length(y));
21.     for i = 1:length(y)
22.         residual(i) = p_y(i) - y(i);
23.     end
24.     fprintf("Residual: ");
25.     disp(residual);
26.     sum = 0;
27.     for i = 1:length(residual)
28.         sum = sum + residual(i) * residual(i);
29.     end
30.     fprintf("Sum of residual squares: %.2f\n", sum);
31.  
32.     figure;
33.     bar(x, residual);
34.     title('Residual values');
35.     xlabel('x');
36.     ylabel('y');
37. end
38.  
39. % Print function
40. function [] = print(p, p_y)    
41.     fprintf("Parameters:\n");
42.     disp(p);
43.     fprintf("Model values:\n ");
44.     disp(p_y);
45. end
46.  
47. % Read the data
48. x = [];
49. y = [];
50. fin = fopen("data.txt", 'r');
51. while ~feof(fin)
52.     line = fgetl(fin);
53.     cells = strsplit(line, ',');
54.     x = [x, str2double(cells{1})];
55.     y = [y, str2double(cells{2})];
56. end
57.  
58. % Curve fitting
59. x_plot = linspace(min(x) - 30, max(x) + 30, 100);
60. x_pred = [12 16 21];
61. y_true = [28.4 28.5 29.5];
62. for j = 1:3
63.     if j == 1
64.         p = polyfit(x, y, 1); % Estimate function parameters
65.         p_y = polyval(p, x); % Estimate function value
66.         y_pred = polyval(p, x_pred);
67.         f = p(1)*x_plot + p(2); % Estimate function
68.         draw(x, x_plot, y, f, p_y, x_pred, y_pred, y_true);
69.         print(p, p_y);
70.         res(x, p_y, y);
71.     elseif j == 2
72.         p = polyfit(x, y, 2); % Estimate function parameters
73.         p_y = polyval(p, x); % Estimate function value
74.         y_pred = polyval(p, x_pred);
75.         f = p(1)*x_plot.^2 + p(2)*x_plot + p(3); % Estimate function
76.         draw(x, x_plot, y, f, p_y, x_pred, y_pred, y_true);
77.         print(p, p_y);
78.         res(x, p_y, y);
79.     else
80.         p = polyfit(x, y, 3); % Estimate function parameters
81.         p_y = polyval(p, x); % Estimate function value
82.         y_pred = polyval(p, x_pred);
83.         f = p(1)*x_plot.^3 + p(2)*x_plot.^2 + p(3)*x_plot + p(4); % Estimate function
84.         draw(x, x_plot, y, f, p_y, x_pred, y_pred, y_true);
85.         print(p, p_y);
86.         res(x, p_y, y);
87.     end
88. end
89.  
90. % Fminsearch
91. % Initial guess
92. x0 = [10, 2];
93.  
94. % The square of sum of residual
95. func = @(params) sum((params(1)*(1 - exp(-params(2)*x)) - y).^2);
96.  
97. % fminsearch implement
98. optimal_params = fminsearch(func, x0);
99.  
100. % Create the estimate function with the result parameter by the fminsearch
101. optimal_func = optimal_params(1)*(1 - exp(-optimal_params(2)*x_plot));
102. optimal_func_pred = optimal_params(1)*(1 - exp(-optimal_params(2)*x_pred));
103.  
104. % Calculate the values with respect to the result parameter
105. p_y = optimal_params(1)*(1 - exp(-optimal_params(2)*x));
106.  
107. %Plot the graph
108. draw(x, x_plot, y, optimal_func, p_y, x_pred, optimal_func_pred, y_true);
109.  
110. %Print out value and the residual
111. print(optimal_params, p_y);
112. res(x, p_y, y);
113.  
114. % Terminate section
115. fclose(fin);