# [R] Problems with integrate

Uwe Ligges ligges at statistik.uni-dortmund.de
Sat Dec 10 18:44:27 CET 2005

```Tolga Uzuner wrote:

> Hi,
>
> Having a weird problem with the integrate function.
>
> I have a function which calculates a loss density: I'd like to integrate
> it to get the distribution.
>
> The loss density function is:
>
> lossdensity<-function(p,Beta,R=0.4){
> # the second derivative of the PDF
> # p is the default probability of the pool at which we are evaluating
> the lossdensity
> # Beta is the correlation with the market factor as a function of K
> # R is the recovery rate
> K=p
> C=qnorm(p) # default threshold for the pool
> A=1/Beta(K)*(C-sqrt(1-Beta(K)^2)*qnorm(K/(1-R)))
> B=qnorm(K/(1-R))
> alpha0=dnorm(A)*sqrt(1-Beta(K)^2)/(dnorm(B)*Beta(K)*(1-R))
> alpha10=-2*dnorm(A)*(C*Beta(K)-A)/(Beta(K)*(1-Beta(K)^2))
> alpha11=(1-R)*dnorm(A)*dnorm(B)*(Beta(K)-A/Beta(K)*(C*Beta(K)-A))/((1-Beta(K)^2)^1.5)
> alpha20=(1-R)*dnorm(A)*dnorm(B)/sqrt(1-Beta(K)^2)
> return(alpha0+(alpha10+alpha11*fd(K,Beta))*fd(K,Beta)+alpha20*fd2(K,Beta))
> }
>
> Beta needs to be passed in as a function:
>
> Beta<-splinefun(c(.03,.06,.09,.12,.22,.6),c(.117,.248,.35,.426,.603,1),method="natural")
>
> Now, when I try out lossdensity, it works fine:
>
>  > lossdensity(.02,Beta,.4)
> [1] 0.1444915
>
> I defined lossdistribution as:
>
> lossdistribution<-function(p,Beta,R=0.4){
> integrate(function(x)lossdensity(x,Beta,R),0.0000000000001,p)}
>
>
> But this gives a weird error:
>
>  > lossdistribution(0.1,Beta)
> Error in "[<-"(`*tmp*`, k, i, value = c(4.29850518211428, 0, 0, 0, 0,  :
>         number of items to replace is not a multiple of replacement length
>

1. We do not have fd and fd2, hence not reproducible.
indentation).
3. What does traceback() tell us?
4. Are you sure lossdensity works on vectors?

Uwe Ligges

> What's interesting is, if I use the area function from library(MASS), it
> works:
>
>  > lossdistribution<-function(p,Beta,R=0.4){
> + area(function(x){lossdensity(x,Beta,R)},0.00001,p)}
>  > lossdistribution(.02,Beta,.4)
> [1] 0.0002177284
>
>
> I could just go ahead and use that, but would like to understand why
> integrate is not working.
>