版塊導(dǎo)航: 正在加載中...

應(yīng)《網(wǎng)絡(luò)安全法》要求，自2017年10月1日起，未進(jìn)行實(shí)名認(rèn)證將不得使用互聯(lián)網(wǎng)跟帖服務(wù)。為保障您的帳號(hào)能夠正常使用，請(qǐng)盡快對(duì)帳號(hào)進(jìn)行手機(jī)號(hào)驗(yàn)證，感謝您的理解與支持！

24小時(shí)熱門版塊排行榜

北京石油化工學(xué)院2026年研究生招生接收調(diào)劑公告

返回列表

我愛信息光學(xué)

新蟲 (初入文壇)

應(yīng)助: 0 (幼兒園)
金幣: 32.5
帖子: 11
在線: 6.5小時(shí)
蟲號(hào): 3075740
注冊(cè): 2014-03-21
專業(yè): 光學(xué)信息獲取與處理

[求助] 最小二乘法編程遇到問(wèn)題求幫助

我用阻尼最小二乘法求解非線性方程組f=acos(d-b)+c最優(yōu)值問(wèn)題，d是自變量，結(jié)果我想得到迭代次數(shù)和a,b,c在迭代中的變化，該怎么修改，求大神幫忙，代碼如下：
clear all;
tic
vx=zeros(1,90);
vy=vx;
vz=vx;
%%初始值
vx(1,

=8;
vy(1,

=8;
vz(1,

=1;
%%加噪聲
vxn=awgn(vx(1,

,0);
vyn=awgn(vy(1,

,0);
vzn=awgn(vz(1,

,0);
%vxn=3*6^(0.5)+randn(90,1);
%vyn=pi/4+randn(90,1);
%vzn=0.5+randn(90,1);
a=((vxn).^2.+(vyn).^2).^(0.5)*3^(0.5)*0.5;
b=atan(vyn./vxn);
c=0.5*vzn;
%%方位角，90個(gè)點(diǎn)
dd=pi/45:pi/45:2*pi;
vr=zeros(1,90);
vro=zeros(1,90);
for i=1:1:90
vr(i)=a(i)*cos(dd(i)-b(i))+c(i);
vro(i)=4*6^(0.5)*cos(dd(i)-pi/4)+0.5;
end
%%
% 計(jì)算函數(shù)f的雅克比矩陣，是解析式
syms a b c y x real;
y=a*cos(x-b)+c;
Jsym=jacobian(y,[a b c])
% 擬合用數(shù)據(jù)。參見《數(shù)學(xué)試驗(yàn)》，p190，例2
data_1=dd;
obs_1=vr;
% 2. LM算法
% 初始
%a0=3*6^(0.5);
%b0=pi/4;
%c0=0.5;
a0=1;
b0=8;
c0=0;
y_init = a0*cos(data_1(1)-b0)+c0;
% 數(shù)據(jù)個(gè)數(shù)
Ndata=length(obs_1);
% 參數(shù)維數(shù)
Nparams=3;
% 迭代最大次數(shù)
n_iters=50;
% LM算法的阻尼系數(shù)初值
lamda=0.1;
% step1: 變量賦值
updateJ=1;
a_est=a0;
b_est=b0;
c_est=c0;
% step2: 迭代
for it=1:n_iters
if updateJ==1
      % 根據(jù)當(dāng)前估計(jì)值，計(jì)算雅克比矩陣
      J=zeros(Ndata,Nparams);
      for i=1:length(data_1)
         %J(i,

=[exp(-b_est*data_1(i)) -a_est*data_1(i)*exp(-b_est*data_1(i))];
J(i,

=[cos(data_1(i)-b_est) a_est*sin(data_1(i)-b_est) 1];
      end
      % 根據(jù)當(dāng)前參數(shù)，得到函數(shù)值
      %y_est = a_est*exp(-b_est*data_1);
      y_est=a_est*cos(data_1-b_est)+c_est;
      % 計(jì)算誤差
      d=obs_1-y_est;
      % 計(jì)算（擬）海塞矩陣
      H=J'*J;
      % 若是第一次迭代，計(jì)算誤差
      if it==1

         e=dot(d,d);
      end
end
% 根據(jù)阻尼系數(shù)lamda混合得到H矩陣
H_lm=H+(lamda*eye(Nparams,Nparams));
% 計(jì)算步長(zhǎng)dp，并根據(jù)步長(zhǎng)計(jì)算新的可能的\參數(shù)估計(jì)值
dp=inv(H_lm)*(J'*d(

);
g = J'*d(

;
a_lm=a_est+dp(1);
b_lm=b_est+dp(2);
c_lm=c_est+dp(3);
% 計(jì)算新的可能估計(jì)值對(duì)應(yīng)的y和計(jì)算殘差e
%y_est_lm = a_lm*exp(-b_lm*data_1);
y_est_lm=a_lm*cos(data_1-b_lm)+c_lm;
d_lm=obs_1-y_est_lm;

e_lm=dot(d_lm,d_lm);
% 根據(jù)誤差，決定如何更新參數(shù)和阻尼系數(shù)
if e_lm<e
      lamda=lamda/10;
      a_est=a_lm;
      b_est=b_lm;
      c_est=c_lm;
      e=e_lm;
      it=it+1
      disp(e);
      updateJ=1;
else
      updateJ=0;
      lamda=lamda*10;

end
end
dis=[a_est/((1+tan(b_est)^2)^(0.5)*3^(0.5)/2) a_est*tan(b_est)/((1+tan(b_est)^2)^(0.5)*3^(0.5)/2) c_est*2]
toc

%%

vrn=a_est*cos(dd-b_est)+c_est;
plot(dd,vr,'r+');
hold on;
plot(dd,vrn);
plot(dd,vro,'g');
hold off;

回復(fù)此樓

» 猜你喜歡

11408軟件工程求調(diào)劑已經(jīng)有3人回復(fù)
一志愿211 初試270分求調(diào)劑已經(jīng)有8人回復(fù)
300求調(diào)劑，材料科學(xué)英一數(shù)二已經(jīng)有6人回復(fù)
學(xué)碩274求調(diào)劑已經(jīng)有9人回復(fù)
本科雙非材料，跨考一志愿華電085801電氣，283求調(diào)劑，任何專業(yè)都可以已經(jīng)有7人回復(fù)
各位老師好，我的一志愿為北京科技大學(xué)085601材料專碩已經(jīng)有3人回復(fù)
317分一志愿南理工材料工程本科湖工大求調(diào)劑已經(jīng)有7人回復(fù)
081200-11408-276學(xué)碩求調(diào)劑已經(jīng)有5人回復(fù)
316求調(diào)劑已經(jīng)有7人回復(fù)
322求調(diào)劑已經(jīng)有7人回復(fù)

1樓 2015-12-21 00:51:28

已閱回復(fù)此樓關(guān)注TA 給TA發(fā)消息送TA紅花 TA的回帖

manmanbobo

新蟲 (初入文壇)

應(yīng)助: 0 (幼兒園)
金幣: 25
紅花: 1
帖子: 37
在線: 12.5小時(shí)
蟲號(hào): 2965461
注冊(cè): 2014-02-09
專業(yè): 地理信息系統(tǒng)

應(yīng)該不難的，你就改個(gè)代碼就行了。

贊一下

回復(fù)此樓

2樓2017-10-11 09:19:21