As the y and n defined within the function, hence it is not requred to consier as input arguments.
function sol= polyreg(x)
year=[1975;1980;1985;1990;1995;2000;2005;2010;2015;2019];
populationjpnworldbank=[111.94;116.78;120.75;123.54;125.44;126.84;127.77;128.07;127.14;126.26];
x=year;
y=populationjpnworldbank;
n=2;
if size(x,2) > size(x,1)
x = x';
end
if size(y,2) > size(y,1)
y = y';
end
nData = size(x,1); % Number of data
nUnks = size(x,2); % Number of unknown
% Number of unknowns initialize
A=zeros(n^nUnks,nUnks);
%% All possible combination of x
for i=1:nUnks
A(:,i)=floor(mod((1:n^nUnks)/n^(i-1),n));
end
A=[A;eye(nUnks)*n];
A=A(sum(A,2)<=n,:);
%% Construct array to solve
legend=cell(size(A,1),1);
X=zeros(size(x,1),size(A,1));
for i=1:size(A,1)
currentPos=find(A(i,:));
currentLegend='';
expression = '';
for j=1:length(currentPos)
if j==1
currentLegend=[currentLegend,'x',num2str(currentPos(j))];
expression=[expression,'x(:,',num2str(currentPos(j)),')'];
if A(i,currentPos(j)) > 1
currentLegend=[currentLegend,'.^',num2str(A(i,currentPos(j)))];
expression=[expression,'.^',num2str(A(i,currentPos(j)))];
end
else
currentLegend=[currentLegend,'.*x',num2str(currentPos(j))];
expression=[expression,'.*x(:,',num2str(currentPos(j)),')'];
if A(i,currentPos(j)) > 1
currentLegend=[currentLegend,'.^',num2str(A(i,currentPos(j)))];
expression=[expression,'.^',num2str(A(i,currentPos(j)))];
end
end
end
if isempty(currentPos)
legend{i,1} = '1';
X(:,i)=ones(size(x,1),1);
else
legend{i,1}=currentLegend;
X(:,i)=eval(expression);
end
end
%% Solution
a=(X'*X);
b=(X'*y);
sol.ab=[a,b];
c = a\b;
sol.legends = legend;
sol.constants = c
%% Function
var='@(';
for i=1:size(A,2)
if i==1
var=[var,'x',num2str(i)];
else
var=[var,',x',num2str(i)];
end
end
var =[var,')'];
fn='';
for i=1:length(legend)
if i ==1
fn = [fn,num2str(c(i)),'.*',legend{i}];
else
fn = [fn,'+',num2str(c(i)),'.*',legend{i}];
end
end
%% Outputs
sol.fn =str2func([var,fn]); %Show function
xs = num2cell(x);
ymodel = zeros(size(xs,1),1);
for i=1:size(xs,1)
ymodel(i) = sol.fn(xs{i,:});
end
%Graph plotting
figure(1)
plot(year, populationjpnworldbank,'*')
xlabel('Year')
ylabel('Number of human population(in Million)')
title('Human Population for Japan in World Bank Data (1975-2019)')
hold on
G = [ ones(length(y),1),x,x.^2];
m = (G'*G)\G'*y;
Dtilde = G*m;
plot(x,Dtilde,'r-')
hold off
%Result
avgy = mean(y);
sol.R2=((sum((y-avgy).^2))-(sum((y-ymodel).^2)))/(sum((y-avgy).^2)); %Show R^2
end
Call the function from command window or another script
>> polyreg(2012)
sol =
struct with fields:
ab: [3×4 double]
legends: {3×1 cell}
constants: [3×1 double]
ans =
struct with fields:
ab: [3×4 double]
legends: {3×1 cell}
constants: [3×1 double]
fn: @(x1)59.0174.*x1+-59118.5762.*1+-0.014697.*x1.^2
R2: 0.9514
