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)
>>[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
>>
>>
>>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.
>>
>>Thanks in advance,
>>Tolga
>>
>>______________________________________________
>>R-help@stat.math.ethz.ch mailing list
>>https://stat.ethz.ch/mailman/listinfo/r-help
>>PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
>>
>>
>>
>
>
>
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)
[1] 0.1444915
 > lossdistribution(.02,Beta,.4)
[1] 0.0002221396
 > lossdistribution(.2,Beta,.4)
[1] 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
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html Received on Sun Dec 11 05:21:10 2005

This archive was generated by hypermail 2.1.8 : Sun 11 Dec 2005 - 14:44:54 EST