It's likely to be quite slow, but it may be possible to take this on numerically. For performance reasons MATLAB insists that integrands be "vectorized", i.e. capable of accepting arrays of inputs and returning arrays of outputs, so we will need to address this. First let me assume that f and p are not coded in a way that would permit me to pass in an array and expect them to return an array of the same size.
gscalar = @(x,y)integral(@(v)arrayfun(@(v)(1 - f(x,y,v)).*p(v),v),y,inf);
g = @(x,y)arrayfun(gscalar,x,y);
T = integral2(@(x,y)exp(g(x,y)),0,inf,@(x)x,inf)
First we defined a gscalar function that accepts scalar x and y and evaluates the innermost integral. Then we used arrayfun to "vectorize" the gscalar function, creating a function that could accept arrays of x and y values and calculates the integrals element-wise. Finally integral2 calculates the outer double integral.
You can tack on 'AbsTol' and 'RelTol' arguments to the integral and integral2 calls.
In case you (or someone reading this) are not familiary with arrayfun, arrayfun is a "loop in a box" sort of thing. Think of
as if it were
z = zeros(size(x));
for k = 1:numel(z)
z(k) = f(x(k),y(k));
but arrayfun should be faster than writing the loop in MATLAB.
You might be able to drop the use of arrayfun in the definition of gscalar if f and p are coded in such a way that they can handle an array of v inputs combined with scalar x and y inputs and deliver the correct value.