How to get the result of int from -inf to inf

I'm a beginner with Matlab, and I'm trying to solve the following problem

h = normpdf(x)./normcdf(y.*x./sqrt(1+2.*y.^2));
E(y) = int(h,x,-inf, inf)

for more, I would like to plot

I(y) = 1+(2.*y.^2./sqrt((1+2.*y.^2).*(2.*pi).^3)).*E(y);

against y. Please help me. Thank you so much.

Answer 1

[Walter Roberson]

syms x E(y) I(y)
h = normpdf(x)./normcdf(y.*x./sqrt(1+2.*y.^2))
E(y) = int(h,x,-inf, inf)
I(y) = 1+(2.*y.^2./sqrt((1+2.*y.^2).*(2.*pi).^3)).*E(y)
fplot(I)

However, in practice the fplot() will be very very very slow. The equations you have given do not have a nice closed form for the expectation or the "I", so every integral has to be computed numerically.

Answer 2

[David Goodmanson]

Hi P^2,

Interesting problem. The integral is a function of y/sqrt(1+2y^2) which I imaginatively called 'a'. For -inf < y < inf, the range of 'a' is (-1/sqrt(2)) < a < (1/sqrt(2)). This is good because I believe the integral only converges for abs(a) < 1. The integral is an even function of 'a' and hence an even function of y.

I don't have the normcdf function but regular Matlab has the error function and

normcdf(x) = (1/2)(1+erf(x/sqrt(2)))

After some algebra, the following code finds E as a function of 'a' by means of convolution, resulting in 29251 array points. The calculation goes up to a = +-.91 even though the definition of a using y stops at +-.707. I compared the results at a few points with straight numerical integration of the initial function using 'integral' and there is agreement within a few parts in 10^7.

After that is done, the approach is to define 'a' in terms of y, and then interpolate the E array as a function of 'a'. For a y array of 20,000 points, the whole process takes about 2 sec on my PC. The second plot shows y vs. the function I. I don't know if this is in line with what you are expecting.

.

% f(a) = A Int{-inf,inf} exp(-x^2/2) / [ A Int{-inf,ax} exp(-u^2/2) du ] dx 
% f(0) = 2 

A = (1/sqrt(2*pi));
% adjustable parameters, result should be fairly insensitive to these
F = 1;  G = 1;
lambda = 1/2;  beta = 1.2;           % need 0<lambda<1 and 1<beta<2
mina = 1e-7;  maxa = .99*(1/sqrt(beta));    % need max(a)<1/sqrt(beta)
N = 100000;

z = linspace(-120/lambda,6-log(beta*F),N);
v = linspace(-120/lambda,6-log(beta*F),N);
delv = v(2)-v(1);
% define   x^2/2 = F*exp(v)   and   1/a^2 - beta = G*exp(w)
R = exp(-F*G*exp(z)).*exp(lambda*z/2);
S1 = exp((1-beta)*F*exp(v)).*exp((1-lambda)*v/2);
S2 = 1/2 + (1/2)*erf(sqrt(F)*exp(v/2));
S3 = (1/2)*erfcx(sqrt(F)*exp(v/2));
S = S1./(S2.*S3);
E = delv*conv(R,fliplr(S));      % correlation
w = delv*(-N+1:N-1);             % lags for correlation
a = 1./sqrt(G*exp(w)+beta);
E = A*sqrt(F/2)*exp(-lambda*w/2).*(1./a).*E;
ind = (a >= mina) & (a <= maxa);
avec = [-a(ind) 0 fliplr(a(ind))];  % symmetric in a
Evec = [ E(ind) 2 fliplr(E(ind))];
% result of main calculation
figure(1);  plot(avec,Evec);  grid on 

% find I vs. y
y = -10:.001:10;
a = y./sqrt(1+2*y.^2);
E = interp1(avec,Evec,a);
I = 1+(2.*y.^2./sqrt((1+2.*y.^2).*(2.*pi).^3)).*E;
figure(2);  plot(y,I);   grid on

% compare results to check
a = .6;
E1 = interp1(avec,Evec,a)
E2 = integral(@(x) efun(x,a), -inf,inf)
E1-E2

function In = efun(x,a)
A = (1/sqrt(2*pi));
In = zeros(size(x));
ind = (x>=0);
In(ind) = A*exp(-(x(ind).^2)/2)./((1/2)*(1+erf((x(ind)*a)/sqrt(2))));
In(~ind) = A*exp(-(1-a^2)*(x(~ind).^2)/2) ...
    ./((1/2)*erfcx((abs(x(~ind))*a)/sqrt(2)));
end