/*to search for maximum in MP problem*/
new;
library optmum,pgraph;

clear d2, X, y;

proc rssc(theta);
	local mmu, msig,Z,r,k, e;
	k = cols(X);
	mmu = theta[1];
	msig = theta[2];
	z = 1./(1+exp(-(d2[.,27]-mmu)/msig));
	X[.,26] = z;
	r = qr(X~y);
/*	{r, e} = qrep(X~y, 
		( 1|zeros(cols(X)-1,1)|(-1) )
					); @  more accurate @*/
	retp( 0.5*rows(X)*ln( r[k+1,k+1]^2 ) );
endp;


load d[]=c:\gausswin\x\BoE\dati1099.csv;
kT=rows(d);
kT "221*44=" 221*43;
d2 = reshape(d,221,43);
"first row";
d2[1,.];
"last row";
d2[221,.];

d2=d2[13:221,.];
y = d2[.,1];
n = rows(d2);
mmu = 0.12; msig=0.21;
z = 1./(1+exp(-(d2[.,27]-mmu)/msig));"mmu " mmu "msig " msig;"z[1] " z[1];
X = ones(n,1)~d2[.,2:25]~z~d2[.,30:43];

output file = exoptmum.out reset;

{theta,f,g,retcode}=optprt(optmum(&rssc, mmu|msig));
mmu = theta[1];
msig = theta[2];
z = 1./(1+exp(-(d2[.,27]-mmu)/msig));
_ptek = "Z.tkf";
title("Z");
xy(seqa(1,1,n),z);

varm = @(2*f/(n*n))*@inv(_opfhess);
se = sqrt(diag(varm));
t = theta./se;
"coeff se t";
theta~se~t;

output off;