How to find maximum from fminunc?

조회 수: 5 (최근 30일)
Boonanan Lechayakittikorn
Boonanan Lechayakittikorn 2018년 1월 30일
편집: Walter Roberson 2018년 1월 30일
Hello everyone, i have problem about to find parameter of beta(0) and beta(1). i have data from generate. My problem is when i input value for generate X and y (value input are stdevu stdevex etc.) when i run code by stdevu = 0.1 and stdevex = 0.1 etc. output likely stdevu = 0.1 and stdevex = 0.3 but i change stdevu = 0.1 and stdevex = 0.1 etc. output unlikely stdevu = 0.3 and stdevex = 0.1 etc. please help me. My code is
function [EX,X,y,W,z,lam,u,uN,vN,EX1] = GEN(XS,stdevu,stdevex,n,p,col)
b = ones(col,1);
uN = normrnd(0,1,[1 n]);
vN = normrnd(0,1,[p n]);
u = stdevu*uN;
%EX1 = normrnd(0,stdevex,n,p);
EX1 = stdevex*vN;
EX = [zeros(n,1) EX1'];
covEX = cov(EX1);
X = XS + EX;
y = (XS*b)+ u';
X_seq=ones(n,1);
for j = 1: col
C=X(:,j) - mean(X(:,j),1);
X_seq=[X_seq C];
end
X_seq(:,1:2)=[];
W = X_seq;
z = y -mean(y);
lam = (stdevex^2)/(stdevu^2);
end .........................................................................................................
clc;
rng(156);
global stdevu stdevex Xdiag V vN s rcovar XS1 y n u uN W z lam ;
stdevu = 0.1;
stdevex = 0.3;
n = 10;
p = 1; %%%%%
col = p+1; %%%%%
it = 1;
tmsedls = 0;
XS1 = normrnd(2,1,n,p);%%%%%
XS = [ones(n,1) XS1];%%%%%
rcovar = diag(stdevex*ones(n*p,1));%%%%%%%%
s = ones(n*p,1);
msedls = 0;
[EX,X,y,W,z,lam,u,uN,vN,EX1] = GEN(XS,stdevu,stdevex,n,p,col);
Xdiag = diag(X(:,2));%%%%%
%V = diag(EX1(:,1));%%%%%%%%%%%%
V = diag(vN);%%%%%%%%%%%%
options = optimoptions(@fminunc,'Algorithm','quasi-newton');
b= [1;1];%beta initial value with b0
f = -(((stdevu*y'*uN'+s'*Xdiag*rcovar*V*s)/n)-(y-b(1)*s-XS1*b(2))'*(y-b(1)*s-XS1*b(2))/(2*n)- (s'*V*rcovar^2*V*s)/(2*n))
[b fval] = fminunc(@(b)-(((stdevu*y'*uN'+s'*Xdiag*rcovar*V*s)/n)-(y-b(1)*s-XS1*b(2))'*(y-b(1)*s- XS1*b(2))/(2*n)-(s'*V*rcovar^2*V*s)/(2*n)),b,options);

답변 (0개)

카테고리

Help CenterFile Exchange에서 Get Started with Curve Fitting Toolbox에 대해 자세히 알아보기

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by