library("metrics")
library("plot")
;
randomize(0)
n = 100
p = 4
beta1 = #(1, 2, 3, 0)
beta1 = beta1/sqrt(beta1'*beta1) ;the first true index
beta2 = #(1,(-2), 1, 1)
beta2 = beta2/sqrt(beta2'*beta2) ; the second true index
x = normal(n, 4) ; independent variable
m =(x*beta1).^2 +(x*beta2).^2 ; true function
e = 0.2*normal(n,1) ; error term
B = beta1~beta2 ; multi-index
h = 1.5*trans(5*n^(-1/(1+p))*sqrt(var(x))) ; bandwidth
D = 2 ;
logi = 0
y = m + e
MI = MIregest(x, y, D, h, logi) ; call of Multi Index regression quantlet
xp = x*MI ; determine projected index
Bhat = MI[,1:D]
esterror = Bhat'*(diag(matrix(p)) - B*B')*Bhat ; compute the accuracy of the projection
esterror
// -- display the multi index model in a D*D plot matrix
prin1 = xp[,1]~y ; the first projection
prin2 = xp[,2]~y ; the second projection
prin12 = xp[,1:2]~y ; the 1-2 projection
dp = createdisplay(D,D)
show(dp, 1, 1, prin1)
setgopt(dp, 1, 1, "title", "First projection")
show(dp, 2, 2, prin2)
setgopt(dp, 2, 2, "title", "Second projection")
show(dp, 2, 1, prin12)
setgopt(dp, 2, 1, "title", "rotate")