Using Quad2d with indicator function

Hi all, I'm having trouble using quad2d (or integral2). I need to find the integral of an empirical copula (see http://en.wikipedia.org/wiki/Copula_%28probability_theory%29#Empirical_copulas) but I find it hard to properly vectorize the function. What I now have is: quad2d(@ecopula,0,1,0,1)

function N= ecopula (X,Y)

load Data_Frees_Valdez.mat
[r2, c2] = size([LOG_LOSS LOG_ALEA]);

C=bivempcop([LOG_LOSS LOG_ALEA]);  

g=zeros(length(X),1);
for index=1:r2
    i = find(X > C(index,1)); 
    for index2 = 1:numel(i)
        if Y(i(index2)) > C(index,2) 
            g(i(index2))=g(i(index2))+1; 
        end
    end
end
result=g./r2; 

N = result;

bivempcop gives me the empirical distribution of the marginal variables and is a [r2,2] double. But when I run this code I get the following error message:

??? Attempted to access g(15); index out of bounds because numel(g)=14.

Error in ==> ecopula at 16 g(i(index2))=g(i(index2))+1; %Indien dit zo is tel er 1 bij (indicator functie)

Error in ==> quad2d>tensor at 355 Z = FUN(X,Y); NFE = NFE + 1;

Error in ==> quad2d at 169 [Qsub,esub] = tensor(thetaL,thetaR,phiB,phiT);

My suspicion is that quad2d passes matrices instead of vectors, but I'm not sure how this exactly works. Any suggestions on how to fix this? thanks

Answer 1

[Mike Hosea]

QUAD2D sends in m-by-m matrices, which you are required to evaluate "elementwise", i.e. if Z = f(X,Y), then Z is the same size as X (and Y), and Z(i,j) = f(X(i,j),Y(i,j)).

The LENGTH function was really only meant for vectors. It is defined on matrices, and even on N-D arrays, as the maximum non-zero dimension length. So, for 14x14 matrix X, length(X) is 14. FIND, on the other hand, linearizes arrays. So an m-by-n input to FIND will result in a vector output of FIND that is between 0 and m*n elements in length. You could do this:

g = zeros(numel(X),1)

But you also need to make sure that the result N has the same size as X before the function returns, so I think what you need is

g = zeros(size(X));

You don't have to change the way you index g in what follows.

Rather than load the data each time, you may want to take advantage of the PERSISTENT feature

persistent r2 c2 C 
if isempty(C)
    load('Data_Frees_Valdez.mat','LOG_LOSS','LOG_ALEA');
    [r2, c2] = size([LOG_LOSS LOG_ALEA]);
    C=bivempcop([LOG_LOSS LOG_ALEA]);
end

This will prevent the function from doing those computations every time it is called. If you change the data set, just type

>> clear ecopula

to force the function to re-read the data on the next call.

Now, moving on from that detail, QUAD2D is not a great method for integrating functions with jump discontinuities. It is really intended for smooth functions with or without edge singularities. I would recommend using the newer INTEGRAL2 function with the 'iterated' method. This integration method handles jump discontinuities much more efficiently. The problem with QUAD2D, and with INTEGRAL2 using the 'tiled' method, is that the code can't really refine the mesh only perpendicularly to the discontinuity. It ends up refining in both x and y directions, which can create a lot of function evaluations. Normally we suggest using the 'tiled' method anyway but after splitting the region so jump discontinuities only occur on integral boundaries, but that's not practical when the function is fundamentally piecewise-defined with many pieces. -- Mike