How to perform recursive integration

Hi,

I am trying to evaluate a function defined by a recursive integral (the math relates to renewal theory. The function is the cumulative probability density function for the ith failure in a renewal process - using a Weibull distribution in the example code).

The recursive series of functions are defined by:

Here is a code example (I want to evaluate F2(0.5) for a weibull distribution with shape = 2 and scale = 1):

wb = @(x) wblpdf(x,2,1);

cdfn(wb,2,0.5)

function y = cdfn(f,i,t)

if i==0

y=1;

else

fun = @(x) f(x)*cdfn(f,i-1,t-x);

y = integral(fun,0,t);

end

end

Edit: I have attached a piece of code that does the job (at least for values of Shape > 1) using 'brute force' integration (= treating it as a discrete sum). The final output is the renewal function:

Answer 1

[Torsten]

syms a b x t positive
syms F0
f(x) = b/a*(x/a)^(b-1)*exp(-(x/a)^b);

f(x) = 

F0(t) = 1;
F1(t) = int(f(x)*F0(t-x),x,0,t)

F1(t) = 

F2(t,x) = f(x)*F1(t-x)

F2(t, x) = 

vpaintegral(subs(F2(t,x),[a b t],[1 0.2 0.5]),x,0,0.5)

ans = 

0.32976

[Dyuman Joshi]

Note - requires Symbolic Math Toolbox.

[Ronnie Vang]

Thanks for the fast response - and also thanks for the clarifying comment. I am new to MatLab and do not know anything about Symbolic Math (I will for sure check it out), so I was actually scratching my head and thinking this was not MatLab code :-)

As I read the code, it is not actually a recursive solution, correct? It is a hard coded version only working for i=2?

In my application, I will be looping through i values (with a condition that F(i,1) is larger than some defined error. The value of i can go to 50 or even higher values.

I made it work using 'brute force' integration (treating it as discrete sums), but I would like to take advantage of the MatLab integration function.

[Torsten]

As I read the code, it is not actually a recursive solution, correct? It is a hard coded version only working for i=2?

You wrote you want F2(0.5). So no need for a recursive solution.

In my application, I will be looping through i values (with a condition that F(i,1) is larger than some defined error. The value of i can go to 50 or even higher values.

I don't understand "F(i,1) is larger than some defined error". You mean "smaller" ?

I made it work using 'brute force' integration (treating it as discrete sums), but I would like to take advantage of the MatLab integration function.

I doubt you can make work a 50-fold convolution integral reliably.

[Ronnie Vang]

Sorry for being unclear (this is my first question in this forum). The call with i=2 was just an example. I want to be able to call it with any value of i.

The 50 fold case may be extreme, it was just to stress that I don't kow the max value of i.

I attached a 'sort of' working code. It is a 'while loop' running while F(i,1) is larger than a defined error value.

[Torsten]

And do you get the result from above (F2(0.5) ~ 0.32976) with your code ? And 1 in the limit as t -> Inf for larger values of i ?

This code should work in the general case, but it will become slow very soon:

a = 1;
b = 0.2;
f = @(x)b/a*(x./a).^(b-1).*exp(-(x/a).^b);
F{1} = @(t) ones(size(t));
for i = 2:3
    F{i} = @(t) integral(@(x)f(x).*F{i-1}(t-x),0,t,'ArrayValued',true);
end
F{3}(0.5)
F{3}(1000)

[Ronnie Vang]

I’m replying from my phone, so unable to edit or run my code.

I’ll test your code tomorrow. Thanks for your support!

[Torsten]

Maybe of interest for you:

https://www.google.com/url?sa=t&rct=j&q=&esrc=s&source=web&cd=&cad=rja&uact=8&ved=2ahUKEwistZG8lLWCAxVRQPEDHShTBo0QFnoECBAQAQ&url=https%3A%2F%2Facademic.oup.com%2Fcomjnl%2Farticle-pdf%2F21%2F3%2F270%2F1061537%2F210270.pdf&usg=AOvVaw0RKaW_czf_luB0xVLcldsV&opi=89978449

Thus a FORTRAN code seems to exist to solve your problem efficiently.

[Ronnie Vang]

And do you get the result from above (F2(0.5) ~ 0.32976) with your code ? And 1 in the limit as t -> Inf for larger values of i ?

Yes, with a shape of 0.2 F2(0.5) = 0.3298. And yes, Fi will approach 1 in the limit of t -> inf for any value of i (Fi is a cumulative distribution function).

Values for i=2 can be done with no recursion (this should be the equivalent to your 'symbolic math example'):

shape = 0.2;
t = 0.5;
fun = @(x) wblpdf(x,1,shape).*wblcdf(t-x,1,shape)

fun = function_handle with value:
    @(x)wblpdf(x,1,shape).*wblcdf(t-x,1,shape)

integral(fun,0,t)

ans = 0.3298

I cannot get your other code to complete. So guess you are right that it will get very slow!

The reference you shared seems to be a bit out of my comfort zone (I am a physicist :-)), so I guess I will stick with the code I have for now.

[Ronnie Vang]

Actually, you provided the answer to my question (I just needed to read what you actually did in your code).

I had not set the ArrayValued property to true. With this my initial example code now works. It is very slow for low values of shape,but that is not a big issue formy use. I am only interested in values of Shape > 1.

Thanks for the support!

wb = @(x) wblpdf(x,1,2);

cdfn(wb,2,0.5)

ans = 0.0094

function y = cdfn(f,i,t)
    if i==0
        y=1;
    else
        fun = @(x) f(x)*cdfn(f,i-1,t-x);
        y = integral(fun,0,t,'ArrayValued',true);
    end
end