/*****************************************************************************
 demo_mata: test & demo for nwspgr, mata version
 integrate function using sparse grids and simulation
*****************************************************************************/

mata:
mata clear
mata set matastrict off

// define g
real colvector g(real matrix x){
	gD = exp(-(x/2):^2/2)/2/sqrt(2*pi() )
	return(exp(rowsum(log(gD))))
}

D = 10    // dimensions
maxk = 4  // max. accuracy level (pol. exactness wil be 2k-1)

// true integral
trueval=(normal(.5)-.5):^D

for (k=2;k<=maxk;k++){
  // integration on sparse grids
  nw = nwspgr("KPU", D, k)
  x = nw[.,1..D]
  w = nw[.,D+1]
  SGappr = g(x)'*w
  SGerror = sqrt((SGappr - trueval):^2):/trueval
  
  // simulation with the same number of nodes, 1000 simulation repetitions
  numnodes = rows(w)
  sim = J(1000,1,.);
  for (r=1;r<=1000;r++){
      x = uniform(numnodes,D)
      sim[r] = mean(g(x))
  }
  Simerror = sqrt(mean((sim:-trueval):^2))/trueval;
  printf("D=%2.0f, k=%2.0f (nodes=%4.0f): SG error=%8.5f, Sim. error=%8.5f\n", D,k,numnodes,SGerror,Simerror )
}

end