Piece wise linear regression from two dimensional data – multiple break points

This Matlab program is commented in Swedish. I give no guarantee that it is working
since it is a long time ago since I wrote it.

Matlab:


% programmet indelar mätpunkter i n över k kombinationer av delmängder
% och returnerar den kombination av korrelationer för linjära
% regressionslinjer, som har den minsta standardavvikelsen.
% Dvs. Programmet försöker hitta en kombination av linjer
% som alla passar lika bra på var sitt avsnitt av mätdata.

lines = 4; % antal linjer
gemensamma_punkter = 0; % 1 = gemensamma punkter beaktas, 0 gemensamma punkter beaktas inte

contactangle = [0.44, 71.1747411262677;
0.46, 70.6802138431153;
0.5, 68.7976902726958;
0.56, 65.9770432272691;
0.64, 63.6548931533462;
0.72, 61.3545600435929;
0.8, 59.0775645191162;
1, 58.6705339485564;
1.5, 56.0490010641405;
2, 54.2248928559115;
3, 53.855355702081;
4, 52.4002296826947;
5, 51.6887751091207;
7, 50.2296697722891;
10, 50.2296697722891;
15, 49.8878105268492;
20, 49.8878105268492;
30, 48.4185494511706;
40, 46.9321517393474;
60, 43.9088037168348;
90, 42.0750220508436;
];

%---------------------------------------------------------------------
x_data = contactangle(:,1);
x = contactangle(:,1);
y = contactangle(:,2);

% räkna medelvärde av punkter med identiska x-värden
Y = [];
X = [];
for row = 1:length(x)
index = find(x==x(row));
if index(1,1) == row
Y(row) = mean(y(index));
X(row) = x(row);
else
Y(row)=NaN;
X(row)=NaN;
end
index = [];
end
contactangle = contactangle(find(not(isnan(Y))),:);
% räkna medelvärde av punkter med identiska x-värden, slut

x = contactangle(:,1);
y = contactangle(:,2);

x = x'; % transponering
y = y'; % transponering

%x = [1, 1.5, 2, 3, 4, 5];% 6 stycken mätpunkter x = tiden
%y = [71.79475502, 69.83880402, 69.83880402, 69.83880402, 68.37224497, 68.37224497];% y = kontakt vinkel

k =[];
k = lines - 1; % k = 3-1 motsvarar 3 linjer
n = length(x)-1; % antal mätdata - 1
cmb =[];
kombination = [];
kombination = nchoosek(1:n,k); %listning av kombinationer
cmb = nchoosek(n,k); %antal kombinationer

index = [1:n];
subsets =[];
partitions = [];
partitions(1:cmb,1:n)=0;
counter =0;
for row = 1:cmb
for column = 1:k
counter =counter +1;
for i=1:n
if kombination(row,column)==index(i)
partitions(row,i)=counter;%kombination(row,column); %partitions(row,i)=kombination(row,column);
end
end
end
counter = 0;
end
partitions; % korrekt

presubsets = [];
onescolumn = [];
onescolumn(1:cmb,1)=1;
presubsets = [onescolumn partitions]; %antalet positioner för partitions strecken är ett mindre än antalet mätdata. Därför lägger man till en column med 1 i början.

for row = 1:cmb
for column = 1:length(x)
if and(column==1,presubsets(row,column)>0)
presubsets(row,column)=1;
end
if and(column>1,presubsets(row,column)>0)
presubsets(row,column)=presubsets(row,column)+1;
else
if column>1
presubsets(row,column)=presubsets(row,column-1);
end
end
end
end
subsets = presubsets; %
% summan av delmängderna är talserie A124932 i
% Online Encyclopedia of integer sequences
% http://www.research.att.com/~njas/sequences/A124932
% sum(sum(subsets));
index = [];
cmb;
chosen_subset_combinations=[];
appearances = [];
counter = 0;

% det här hade kanske gått att göra enklare från början
for row = 1:cmb
for line=1:lines
index = find(subsets(row,:)==line);
appearances = [appearances length(index)];
index = [];
end
if min(appearances)>1
counter = counter +1;
chosen_subset_combinations(counter,:)=subsets(row,:); % väljer de kombinationer som vars subset innehåller 2 eller flera punkter
end
appearances = [];
end

chosen_subset_combinations; % valda kombinationer av subsets

subsets = [];
subsets = chosen_subset_combinations;

[m,n]=size(subsets);
index=[];
X = [];
Y = [];
R2 = [];
R = [];
r = [];
%R2(m,n) = 0;
for row = 1:m
for column = 1:lines
index = find(subsets(row,:)==column);
if gemensamma_punkter == 1
% beaktande av gemensamma punkter för subsets, början
first_point = min(index);
last_point = max(index);
if and(first_point==1, last_point1, last_point1, last_point==length(x))
index = [(first_point-1) index];
end
end
end
% beaktande av gemensamma punkter för subsets, slut
end
X = x(index);
Y = y(index);
% p = polyfit(X,Y,1);
% Y_estimate = p(1)*X + p(2);
if std(Y)==0 % om mätdata innehåller konstanta platåer(=linjer) skall r = 1
r=1;
else
R = corrcoef(X,Y);
r = R(1,2);
detect = [];
detect = isnan(r); % Det här vilkoret är kanske litet riskabelt och kanske även onödigt.
if detect == 1
r = 1; % corrcoeff NaN värdena och varningsmeddelandena torde härstamma från de
% konstanta segmentena i mätdata. Detta borde utredas.
% Åtgärdat ovan med vilkoret std(Y)==0
end
end
R2(row,column) = r^2;
end
end
R2;

% summera R2 värdena och se vilken kombination som har det högsta sammanlagda R2 värdet
% alternativt se vilken kombination som har lägsta standard avvikelsen.
% Jag vet inte vilket som är bäst.

s=[]; % målfunktionen s
for row=1:m
% s(row) = var(R2(row,1:lines)); % variansen av linjernas korrelation
s(row) = std(R2(row,1:lines)); % standardavvikelsen av linjernas korrelation
% s(row) = sum(R2(row,1:lines)); % summan av linjernas korrelation
end
s';

[C,I] = min(s); % minimering
%[C,I] = max(s); % maximering

combination_with_lowest_standard_deviation_of_R2_values = chosen_subset_combinations(I,:);
combination_with_lowest_standard_deviation_of_R2_values = combination_with_lowest_standard_deviation_of_R2_values'
R2_values_for_chosen_combination = R2(I,:)
% average_of_R2_values_for_chosen_combination=mean(R2(I,:));

Advertisements
This entry was posted in Uncategorized. Bookmark the permalink.