如何准确拟合旋转双曲线?
How to fit rotated hyperbola accurately?
我的实验数据点看起来像一条条双曲线。下面我提供了一个代码(Matlab),它生成“虚拟”数据,与原始数据非常相似:
function [x_out,y_out,alpha1,alpha2,ecK,offsetX,offsetY,branchDirection] = dummyGenerator(mu_alpha)
alpha_range=0.1;
numberPoint2Return=100; % number of points to return
ecK=10.^((rand(1)-0.5)*2*2); % eccentricity-related parameter
% slope of the first asimptote (radians)
alpha1 = ((rand(1)-0.5)*alpha_range+mu_alpha);
% slope of the first asimptote (radians)
alpha2 = -((rand(1)-0.5)*alpha_range+mu_alpha);
beta = pi-abs(alpha1-alpha2); % angle between asimptotes (radians)
branchDirection = datasample([0,1],1); % where branch directed
% up: branchDirection==0;
% down: branchDirection==1;
% generate branch
x = logspace(-3,2,numberPoint2Return*100)'; %over sampling
y = (tan(pi/2-beta)*x+ecK./x);
% rotate branch using branchDirection
theta = -(pi/2-alpha1)-pi*branchDirection;
% get rotation matrix
rotM = [ cos(theta), -sin(theta);
sin(theta), cos(theta) ];
% get rotated coordinates
XY1=[x,y]*rotM;
x1=XY1(:,1); y1=XY1(:,2);
% remove possible Inf
x1(~isfinite(y1))=[];
y1(~isfinite(y1))=[];
% add noise
y1=((rand(numel(y1),1)-0.5)+y1);
% downsampling
%x_out=linspace(min(x1),max(x1),numberPoint2Return)';
x_out=linspace(-10,10,numberPoint2Return)';
y_out=interp1(x1,y1,x_out,'nearest');
% randomize offset
offsetX=(rand(1)-0.5)*50;
offsetY=(rand(1)-0.5)*50;
x_out=x_out+offsetX;
y_out=y_out+offsetY;
end
典型结果如图所示:
数据有以下重要 属性:两条渐近线的斜率来自相同的分布(只是符号不同),所以对于我的拟合,我对 mu_alpha
.[=16 的估计相当不错=]
现在开始有问题的部分。我尝试拟合这些数据点。我的方法的主要思想是找到一个旋转以获得 y=k*x+b/x
形状,然后正好适合它。
我使用以下代码:
function Rsquare = fitFunction(x,y,alpha1,alpha2,ecK,offsetX,offsetY)
R=[];
for branchDirection=[0 1]
% translate back
xt=x-offsetX;
yt=y-offsetY;
% rotate back
theta = (pi/2-alpha1)+pi*branchDirection;
rotM = [ cos(theta), -sin(theta);
sin(theta), cos(theta) ];
XY1=[xt,yt]*rotM;
x1=XY1(:,1); y1=XY1(:,2);
% get fitted values
beta = pi-abs(alpha1-alpha2);
%xf = logspace(-3,2,10^3)';
y1=y1(x1>0);
x1=x1(x1>0);
%x1=x1-min(x1);
xf=sort(x1);
yf=(tan(pi/2-beta)*xf+ecK./xf);
R(end+1)=sum((xf-x1).^2+(yf-y1).^2);
end
Rsquare=min(R);
end
不幸的是,这段代码效果不佳,我经常得到不好的结果,即使我使用已知的(来自模拟的)初始参数也是如此。
你能帮我找到一个好的解决方案吗?
更新:
我找到了解决方案(参见答案),但是
我还有一个小问题——我对a
参数的估计是错误的,有时因为这个原因我没有很好的拟合。
你能提出一些如何从实验点估计 a 的想法吗?
我找到了主要问题(通常是我的大脑)!我不知道双曲线的一般方程。所以我的双曲线方程是:
((x-x0)/a).^2-((y-y0)/b).^2=-1
所以,我们可能不关心符号,那么我可能会使用以下代码:
mu_alpha=pi/6;
[x,y,alpha1,alpha2,ecK,offsetX,offsetY,branchDirection] = dummyGenerator(mu_alpha);
% hyperb=@(alpha,a,x0,y0) tan(alpha)*a*sqrt(((x-x0)/a).^2+1)+y0;
hyperb=@(x,P) tan(P(1))*P(2)*sqrt(((x-P(3))./P(2)).^2+1)+P(4);
cost =@(P) fitFunction(x,y,P);
x0=mean(x);
y0=mean(y);
a=(max(x)-min(x))./20;
P0=[mu_alpha,a,x0,y0];
[P,fval] = fminsearch(cost,P0);
hold all
plot(x,y,'-o')
plot(x,hyperb(x,P))
function Rsquare = fitFunction(x,y,P)
%x=sort(x);
yf=tan(P(1))*P(2)*sqrt(((x-P(3))./P(2)).^2+1)+P(4);
Rsquare=sum((yf-y).^2);
end
P.S。 LaTex 标签对我不起作用
我的实验数据点看起来像一条条双曲线。下面我提供了一个代码(Matlab),它生成“虚拟”数据,与原始数据非常相似:
function [x_out,y_out,alpha1,alpha2,ecK,offsetX,offsetY,branchDirection] = dummyGenerator(mu_alpha)
alpha_range=0.1;
numberPoint2Return=100; % number of points to return
ecK=10.^((rand(1)-0.5)*2*2); % eccentricity-related parameter
% slope of the first asimptote (radians)
alpha1 = ((rand(1)-0.5)*alpha_range+mu_alpha);
% slope of the first asimptote (radians)
alpha2 = -((rand(1)-0.5)*alpha_range+mu_alpha);
beta = pi-abs(alpha1-alpha2); % angle between asimptotes (radians)
branchDirection = datasample([0,1],1); % where branch directed
% up: branchDirection==0;
% down: branchDirection==1;
% generate branch
x = logspace(-3,2,numberPoint2Return*100)'; %over sampling
y = (tan(pi/2-beta)*x+ecK./x);
% rotate branch using branchDirection
theta = -(pi/2-alpha1)-pi*branchDirection;
% get rotation matrix
rotM = [ cos(theta), -sin(theta);
sin(theta), cos(theta) ];
% get rotated coordinates
XY1=[x,y]*rotM;
x1=XY1(:,1); y1=XY1(:,2);
% remove possible Inf
x1(~isfinite(y1))=[];
y1(~isfinite(y1))=[];
% add noise
y1=((rand(numel(y1),1)-0.5)+y1);
% downsampling
%x_out=linspace(min(x1),max(x1),numberPoint2Return)';
x_out=linspace(-10,10,numberPoint2Return)';
y_out=interp1(x1,y1,x_out,'nearest');
% randomize offset
offsetX=(rand(1)-0.5)*50;
offsetY=(rand(1)-0.5)*50;
x_out=x_out+offsetX;
y_out=y_out+offsetY;
end
典型结果如图所示:
数据有以下重要 属性:两条渐近线的斜率来自相同的分布(只是符号不同),所以对于我的拟合,我对 mu_alpha
.[=16 的估计相当不错=]
现在开始有问题的部分。我尝试拟合这些数据点。我的方法的主要思想是找到一个旋转以获得 y=k*x+b/x
形状,然后正好适合它。
我使用以下代码:
function Rsquare = fitFunction(x,y,alpha1,alpha2,ecK,offsetX,offsetY)
R=[];
for branchDirection=[0 1]
% translate back
xt=x-offsetX;
yt=y-offsetY;
% rotate back
theta = (pi/2-alpha1)+pi*branchDirection;
rotM = [ cos(theta), -sin(theta);
sin(theta), cos(theta) ];
XY1=[xt,yt]*rotM;
x1=XY1(:,1); y1=XY1(:,2);
% get fitted values
beta = pi-abs(alpha1-alpha2);
%xf = logspace(-3,2,10^3)';
y1=y1(x1>0);
x1=x1(x1>0);
%x1=x1-min(x1);
xf=sort(x1);
yf=(tan(pi/2-beta)*xf+ecK./xf);
R(end+1)=sum((xf-x1).^2+(yf-y1).^2);
end
Rsquare=min(R);
end
不幸的是,这段代码效果不佳,我经常得到不好的结果,即使我使用已知的(来自模拟的)初始参数也是如此。
你能帮我找到一个好的解决方案吗?
更新:
我找到了解决方案(参见答案),但是
我还有一个小问题——我对a
参数的估计是错误的,有时因为这个原因我没有很好的拟合。
你能提出一些如何从实验点估计 a 的想法吗?
我找到了主要问题(通常是我的大脑)!我不知道双曲线的一般方程。所以我的双曲线方程是:
((x-x0)/a).^2-((y-y0)/b).^2=-1
所以,我们可能不关心符号,那么我可能会使用以下代码:
mu_alpha=pi/6;
[x,y,alpha1,alpha2,ecK,offsetX,offsetY,branchDirection] = dummyGenerator(mu_alpha);
% hyperb=@(alpha,a,x0,y0) tan(alpha)*a*sqrt(((x-x0)/a).^2+1)+y0;
hyperb=@(x,P) tan(P(1))*P(2)*sqrt(((x-P(3))./P(2)).^2+1)+P(4);
cost =@(P) fitFunction(x,y,P);
x0=mean(x);
y0=mean(y);
a=(max(x)-min(x))./20;
P0=[mu_alpha,a,x0,y0];
[P,fval] = fminsearch(cost,P0);
hold all
plot(x,y,'-o')
plot(x,hyperb(x,P))
function Rsquare = fitFunction(x,y,P)
%x=sort(x);
yf=tan(P(1))*P(2)*sqrt(((x-P(3))./P(2)).^2+1)+P(4);
Rsquare=sum((yf-y).^2);
end
P.S。 LaTex 标签对我不起作用