Problem using trapz for obtaining probabilities

I use numerical integration to calculate the the integral as seen below. Because this represents a loss distribution, values for P_D should be in the interval [0,1]. The problem is that e.g. P_D(d=0)>10!

How can I fix this problem?

PD  = .01;
N   = 50;
r   = .1;

    for d=0:N

        f = -10:.01:10;
        CPD = normcdf(norminv(PD,0,1),sqrt(r)*f,sqrt(1-r));
        l = CPD.^(d).* (1-CPD).^(N-d);
        P_D(1,d+1) = nchoosek(N,d) * trapz(f,l);

    end

Answer 1

[Tom Lane]

I don't quite understand what you are trying to do. I suggest you consider d=0 and step through your code to see why the first number it calculates is greater than 1.

I am suspicious that

nchoosek(N,d) * CPD.^(d).* (1-CPD).^(N-d);

has the form of a binomial distribution probability function, but you are integrating it over f (CPD is a function of f) in a way I don't understand.

In the first iteration with d=0, try plot(f,l) right before the call to trapz. You are integrating a thing that looks like a cdf rather than a pdf.