一次式によるベースライン補正について

グラフで得られているx軸0-36までの平均値は変えずに、

最小二乗近似法で得られたオレンジ点線の直線（y = -1.97x+116.08）と「同じ傾きで右肩上がりの直線」を求め、その直線に従いy軸を変化させるにはどのようにしたらいいでしょうか？

もしくは、

x軸が大きくなるにつれて、y軸が高いところを低く、低いところを高く、平均値は変えず全体的に右肩上がりに補正できる方法はないでしょうか？

スペクトルではないので、ベースライン補正なのか分かりませんが、簡易的な一次式で補正できれば考えています。

Answer 1

[Hiro Yoshino]

面積の補正が必用ということで、いただいたコードを少し変えてみました。

利便性を考慮して最初のデータでも動作するようにしました。y が before, y2 が after です。

% Original data
x = (1:36)';
y = [0;134.7407;205.9038;210.5938;148.8053;101.2817;94.7168;89.2617;92.2141;100.3457;103.1461;98.5113;72.87;72.7713;71.103;60.4682;53.5974;58.5912;51.7347;51.1305;50.0097;44.1407;49.0088;45.1054;52.8028;49.5283;68.3816;58.6178;75.5166;67.1995;82.3074;84.7901;103.3022;165.2237;0;0];

plot(x,y);
title("元のデータ");

トレンド除去

idx = find(y>0); % y>0 のインデックス
idx0 = find(y<=0); % y<=0 のインデックス
d = detrend(y(idx),1); % y>0 だけに限定してトレンド除去後を計算

% 全てのデータポイントでトレンド除去を考慮
y2 = zeros(size(y));
y2(idx) = d;
y2(idx0) = y(idx0); % y<=0 の部分には元のデータを保持

トレンド除去後の面積補正

% 面積の補正量
sComp = mean(y(idx)-y2(idx))

sComp = 86.9007

y2(idx) = y2(idx) + sComp; % 面積補正後

% 比較
figure;
tiledlayout(2,2);

nexttile
scatter(x,y,'b');
lsline;
title("Before と LS線")

nexttile
scatter(x,y2,'m');
lsline
title("After と LS線, '0'のためLS線が下がる");

nexttile
plot(x,y,x,y2);
title("2つの比較");

nexttile
plot(x,y-y2);
title(sprintf("面積が等しい証明: mean(y-y2)=%1.3f",mean(y-y2)));
legend("y-y2");

[KT]

返信が遅れて大変申し訳ありません。

トレンド除去の妥当性や面積補正の方法などを考えておりました。

非常に有用なコード作成までご教示ありがとうございました。

Answer 2

[Hiro Yoshino]

回転角を考えて、全体のデータを回転させると良いかと思います。

元のデータ

x = 0:0.1:36;
y = -1.97*x + 116.08;
plot(x,y,'b:');

回転させます

theta = 2 * (pi/2 - atan(-1.97)); % 回転角を計算
R = [cos(theta) -sin(theta); sin(theta) cos(theta)]; % 回転行列
center = repmat([mean(x); mean(y)], 1, length(x)); % 中心座標
s = [x;y] - center; % 原点に移動
s0 = R*s; % 回転
x0y0 = s0 + center; % 元に戻す

plot(x,y,'b:',x0y0(1,:),x0y0(2,:),'r:');

[KT]

非常に分かりやすいご回答ありがとうございます．

このとき，添付した青線の面積（ｙ軸の値の合計）もしくは平均（面積からプロット数36で割った値）は変化しないでしょうか？

最小二乗法では誤差が最小になるように求めるので，もともとの青線の面積と近い値にはなるかもしれませんが，いかがでしょうか．

[Hiro Yoshino]

正直、ご質問がよく理解ができていません。直線の平均値を変えずに傾きを変えることは分かったのですが..何がやりたいか、どこが分からないのか簡潔に教えていただけますか?

最小二乗法は二乗誤差を平均的に最小にする手法です。データによっては面積が一致するでしょうし、全く異なるものになる場合があります。一般化は出来ないかと。

[KT]

青線のデータをx軸が増加するにつれてy軸が増加するように補正する方法が知りたいです。このような青線をベースライン補正することができますか？

まずどのような直線かどうか判別するために近似直線を描きました（必ず右肩下がりの直線になるはずです）。

次に、近似直線を右肩上がりの直線に変換し、その直線の通りにy軸を変換することで補正ができるのではと考えました。

ただし、前提として、青線の面積（ｙ軸の値の合計）は不変とする　というものです。

そのため、最小二乗法から得られた直線では近しい面積にはなるかもしれないが、必ず等しくなるわけではない気がしました。

分かりにくく申し訳ありません。

[Hiro Yoshino]

まず、段階的に確認させてください。

1. 「右肩上がりの直線に変換し」 - 変換する方法 (回転) を示しました。この部分はこの方法で大丈夫ですか？
2. 「その直線の通りに y軸を変換する」 - "直線の通りに" とは数学的にどういう意味ですか?、y軸を変換というのは、どう変換したいのですか？写像ですか、回転ですが、この部分で何がしたいのか具体的に教えてください。
3. 「前提として、青色の面積 (ｙ軸の値の合計) 」 - ここももう少し数学的に記述いただけますか? 積分値を保ったまま、「回転？写像？」。「これを制約条件として xxx という最適化問題を解きたい」位まで落とし込めると、説明がかなり簡単になります。

可能でしたら、やりたいことを実現するためのアルゴリズムを指し示していただけると、MATLAB で実装するにはどうするのか？という議論だけに集中できるので簡単だと思います。

[Hiro Yoshino]

今のベースで少し修正するとこんな感じになります。ご評価あればいただけると、ここからスタートして変更することもできるかと。

% データ作成
x = (0:36)';
y = 2*randn(size(x)) + 10*sin(x*2*pi/20) - 1.97*x + 116.08;
plot(x,y,'o');

[lfit, gof] = fit(x,y,'poly1') % 一次関数をデータにフィット

lfit = 
     Linear model Poly1:
     lfit(x) = p1*x + p2
     Coefficients (with 95% confidence bounds):
       p1 =      -2.142  (-2.358, -1.926)
       p2 =       119.8  (115.2, 124.3)

gof = struct with fields:
           sse: 1.6722e+03
       rsquare: 0.9205
           dfe: 35
    adjrsquare: 0.9182
          rmse: 6.9122

plot(lfit,x,y);

coeffs = coeffvalues(lfit)  % モデルの係数を取得

coeffs = 1×2

-2.1421 119.7550

%theta = 2 * (pi/2 - atan(coeffs(1))) % 回転角を計算
theta = -atan(coeffs(1)) % 回転角を計算

theta = 1.1340

R = [cos(theta) -sin(theta); sin(theta) cos(theta)]; % 回転行列を計算
center = repmat([mean(x); mean(y)], 1, length(x)); % データの中心座標
s = [x,y]' - center; % 原点に移動
s0 = R*s; % 回転
x0y0 = s0 + center; % 元に戻す

plot(x,y,'bo-',x0y0(1,:),x0y0(2,:),'ro-'); % プロット

[Hiro Yoshino]

ベースライン補正 で何となくどんなことをされたいのかは分かりました。

(校正のことなのかなと理解しました)

その上で、ベースライン補正で利用されるアルゴリズムなどあれば、ご紹介いただけると取り掛かり易いです。

[KT]

数学的に議論したいのですが、知識不足で他文献もないため、分かりにくく大変申し訳ありません。

yデータに対し一次式でのベースライン補正を考えたときに、detrend(x,1)が使えるのではないかと考えました。

% x, yデータ
x = (1:36)';
y = [0;134.7407;205.9038;210.5938;148.8053;101.2817;94.7168;89.2617;92.2141;100.3457;103.1461;98.5113;72.87;72.7713;71.103;60.4682;53.5974;58.5912;51.7347;51.1305;50.0097;44.1407;49.0088;45.1054;52.8028;49.5283;68.3816;58.6178;75.5166;67.1995;82.3074;84.7901;103.3022;165.2237;0;0]

y = 36×1

0 134.7407 205.9038 210.5938 148.8053 101.2817 94.7168 89.2617 92.2141 100.3457

% 線形近似　データの傾向として右肩下がりなのか確認のため
X = [ones(length(x),1) x ];
b = X\y

b = 2×1

116.0829 -1.9689

yCalc = X*b;
plot(x,y,'-')
hold on 
plot(x,yCalc,'--')


% 線形トレンド
D = detrend(y,1);
plot(x,D)
plot(x,y-D,":k")
legend("Input Data","Detrended Data","Trend")

y値は正のため、0を含めてしまうとトレンド除去の際にy値0が負に補正されてしまうため、以下では0は除去することを考えました。

x = (1:33)';
y = [134.7407;205.9038;210.5938;148.8053;101.2817;94.7168;89.2617;92.2141;100.3457;103.1461;98.5113;72.87;72.7713;71.103;60.4682;53.5974;58.5912;51.7347;51.1305;50.0097;44.1407;49.0088;45.1054;52.8028;49.5283;68.3816;58.6178;75.5166;67.1995;82.3074;84.7901;103.3022;165.2237]

y = 33×1

134.7407 205.9038 210.5938 148.8053 101.2817 94.7168 89.2617 92.2141 100.3457 103.1461

% 線形近似　データの傾向として右肩下がりなのか確認のため
X = [ones(length(x),1) x ];
b = X\y

b = 2×1

122.2141 -2.0773

yCalc = X*b;
plot(x,y,'-')
hold on 
plot(x,yCalc,'--')

% 線形トレンド
D = detrend(y,1);
plot(x,D)
plot(x,y-D,":k")
legend("Input Data","Detrended Data","Trend")

ただし、これでは積分値が変化することになります。これを変化させないようにトレンド除去する方法はありますでしょうか？

＞「右肩上がりの直線に変換し」 - 変換する方法 (回転) を示しました。この部分はこの方法で大丈夫ですか？

よく分かりました。回転を利用することで、近似曲線を右肩上がりに変換できるので納得しました。ただ、変換した近似曲線にしたがってy値が補正できればと考えています。