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
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. 2. Your code is almost unreadable, please try to help the helpers and make your code readable before posting (by using blanks/spaces and 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. > > Thanks in advance, > Tolga > > ______________________________________________ > R-help at 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
integrate returns a list, not just the value. Try integrate(...)$value See ?integrate for details. On 12/10/05, Tolga Uzuner <tolga at 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 at 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 >