How to find maximum from fminunc

1 vue (au cours des 30 derniers jours)
Boonanan Lechayakittikorn
Boonanan Lechayakittikorn le 29 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);

Réponses (0)

Catégories

En savoir plus sur Get Started with Curve Fitting Toolbox dans Help Center et File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by