# Re: [R] Problems with integrate

From: Tolga Uzuner <tolga_at_coubros.com>
Date: Sun 11 Dec 2005 - 05:11:01 EST

Gabor Grothendieck wrote:

>integrate returns a list, not just the value. Try integrate(...)\$value
>See ?integrate for details.
>
>On 12/10/05, Tolga Uzuner <tolga@coubros.com> 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)
>> 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
>>
>>
>>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)
>> 0.0002177284
>>
>>
>>I could just go ahead and use that, but would like to understand why
>>integrate is not working.
>>
>>Tolga
>>
>>______________________________________________
>>R-help@stat.math.ethz.ch mailing list
>>https://stat.ethz.ch/mailman/listinfo/r-help
>>
>>
>>
>
>
>
Thanks to everyone, vectorizing worked:

lossdensity<-function(ProbVect,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
sapply(ProbVect, function (p) {

```    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))

}
)
}

> lossdensity(.02,Beta,.4)
 0.1444915
> lossdistribution(.02,Beta,.4)
 0.0002221396
> lossdistribution(.2,Beta,.4)
 0.2564789
> lossdistribution
function(p,Beta,R=0.4){
integrate(function(x)lossdensity(x,Beta,R),0.0000000000001,p)\$value}

R-help@stat.math.ethz.ch mailing list