How to find maximum from fminunc?

6 views (last 30 days)
Boonanan Lechayakittikorn
Boonanan Lechayakittikorn on 30 Jan 2018
Edited: Walter Roberson on 30 Jan 2018
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);

Answers (0)

Categories

Find more on Get Started with Curve Fitting Toolbox in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!