0

I am using numerical integration in MATLAB, with one varibale to integrate over but the function also contains a variable number of terms depending on the dimension of my data. Right now this looks like the following for the 2-dimensional case:

for t = 1:T  
   fxt = @(u)  exp(-0.5*(x(t,1)-theta*norminv(u,0,1)).^2) .* ...
         exp(-0.5*(x(t,2) -theta*norminv(u,0,1)).^2);
   f(t) = integral(fxt,1e-4,1-1e-4,'AbsTol',1e-3); 
end

I would like to have this function flexible in the sense that there could be any number of data points in, each in the following term:

exp(-0.5*(x(t,i) -theta*norminv(u,0,1)).^2);

I hope this is understandable.

InfiniteVariance
  • 81
  • 1
  • 8

2 Answers2

0

If x and u have a valid dimension match (vector-vector or array-scalar) for the subtraction, you can put the whole matrix x into the handle and pass it to the integral function using the name-parameter pair ('ArrayValued',true):

fxt = @(u)  exp(-0.5*(x - theta*norminv(u,0,1)).^2) .* ...
            exp(-0.5*(x - theta*norminv(u,0,1)).^2);
f   = integral(fxt,1e-4,1-1e-4,'AbsTol',1e-3,'ArrayValued',true);

[Documentation]

You may need a loop if integral ever passes a vector u into the handle. But in looking at how the integral function is written, the integration nodes are entered as scalars for array-valued functions, so the loop shouldn't be necessary unless some weird dimension-mismatch error is thrown.


Array-Valued Output

In response to the comments below, you could try this function handle:

fx = @(u,t,k)  prod(exp(-0.5*(x(t,1:k)-theta*norminv(u,0,1)).^2),2);

Then your current loop would look like

fx = @(u,t,k) prod(exp(-0.5*(x(t,1:k)-theta*norminv(u,0,1)).^2),2);
k  = 2;
for t = 1:T
  f(t) = integral(@(u)fx(u,t,k),1e-4,1-1e-4,'AbsTol',1e-3,'ArrayValued',true);
end

The ArrayValued flag is needed since x and u will have a dimension mismatch. In this form, another loop would be needed to sweep through the k indexes. However, we can improve this function by skipping the loop altogether since each iterate of the loop is independent by using the ArrayValued mode:

fx = @(u,k) prod(exp(-0.5*(x(:,1:k)-theta*norminv(u,0,1)).^2),2);
k  = 2;
f  = integral(@(u)fx(u,k),1e-4,1-1e-4,'AbsTol',1e-3,'ArrayValued',true);


Vector-Valued Output

If ArrayValued is not desired, which may be the case if the integration requires a lot of subdivisions and a vector-valued u is preferable, you can also try a recursive version of the handle using cell arrays:

% x has size [T,K]
fx = cell(K,1);
fx{1} = @(u,t) exp(-0.5*(x(t,1) - theta*norminv(u,0,1)).^2);
for k = 2:K
    fx{k} = @(u,t) fx{k-1}(u,t).*exp(-0.5*(x(t,k) - theta*norminv(u,0,1)).^2);
end

f(T) = 0;
k    = 2;
for t = 1:T
    f(t) = integral(@(u)fx{k}(u,t),1e-4,1-1e-4,'AbsTol',1e-3);
end
TroyHaskin
  • 8,361
  • 3
  • 22
  • 22
  • Thanks for your effort, I just realized that my code above is wrong, The first x(t,2) must be changed to x(t,1) - so basically the funtion is the product of different terms of various number, each including a different x. – InfiniteVariance Oct 23 '14 at 21:28
  • A different `x` or just a different index? For example, will it always be x(t,i) and x(t,i-1)? – TroyHaskin Oct 23 '14 at 21:30
  • x is Matrix of Dimension T*K and for each t/row I would like to have the function involving all columns, so my above example would be for K = 2, but x could have any dimension. For K = 3 I would like to end up with: exp(-0.5*(x(t,1) - theta*norminv(u,0,1)).^2) .* exp(-0.5*(x(t,2) - theta*norminv(u,0,1)).^2)*exp(-0.5*(x(t,3) - theta*norminv(u,0,1)).^2) .*; – InfiniteVariance Oct 23 '14 at 21:37
0

ThanksTroy but now I run into the follwing:

x = [0.3,0.8;1.5,-0.7];
T = size(x,1);
k  = size(x,2);
theta= 1;
fx = @(u,t,k) prod(exp(-0.5*(x(t,1:k) - theta*norminv(u,0,1))^2));
for t = 1,T
  f(t) = integral(@(u)fx(u,t,k),1e-4,1-1e-4,'AbsTol',1e-3);
end

Error using - Matrix dimensions must agree.

Error in @(u,t,k)prod(exp(-0.5*(x(t,1:k)-theta*norminv(u,0,1))^2))

Error in @(u)fx(u,t,k)

Error in integralCalc/iterateScalarValued (line 314) fx = FUN(t);

Error in integralCalc/vadapt (line 133) [q,errbnd] = iterateScalarValued(u,tinterval,pathlen);

Error in integralCalc (line 76) [q,errbnd] = vadapt(@AtoBInvTransform,interval);

Error in integral (line 89) Q = integralCalc(fun,a,b,opstruct);

InfiniteVariance
  • 81
  • 1
  • 8