function[I, rho, I1, I2, I3] = GaussRomberg(fun, mu, sig, N);
% call with [I, ... ] =GaussRomberg(@fun, ...)
% or [I,...] = GaussRomberg(@(x) abs(x).^alpha, mu, sig, N) e.g.
% See Parameterizing functions or Function functions or Passing a function
%    as an argument or Anonymous functions.
% integrates fun(x)n(x)dx, where n(x) ~ N(mu, sig^2)
p = ((1:N)-.5)/N;
I1 = sum(fun(icdf('normal',p, mu, sig)))/N;
% figure
% plot(p,fun(icdf('normal',p, mu, sig)))
p = ((1:2*N)-.5)/(2*N);
I2 = sum(fun(icdf('normal',p, mu, sig)))/(2*N);
p = ((1:4*N)-.5)/(4*N);
I3 = sum(fun(icdf('normal',p, mu, sig)))/(4*N);
D2 = I2 - I1;
D3 = I3 - I2;
rho = D3/D2;
rho = min(rho, .5);
rho = max(rho, 0);
I = I2 + D3/(1-rho);
return
end
