[Yuima-commits] r688 - pkg/yuima/R

[noreply at r-forge.r-project.org]

Author: yumauehara
Date: 2019-03-04 04:54:09 +0100 (Mon, 04 Mar 2019)
New Revision: 688

Modified:
   pkg/yuima/R/rng.R
Log:
fixed a bug in dbgamma function

Modified: pkg/yuima/R/rng.R
===================================================================
--- pkg/yuima/R/rng.R	2019-02-03 14:49:28 UTC (rev 687)
+++ pkg/yuima/R/rng.R	2019-03-04 03:54:09 UTC (rev 688)
@@ -27,21 +27,25 @@
   if(delta.plus<=0||gamma.plus<=0||delta.minus<=0||gamma.minus<=0)
   stop("All of the parameters are positive.")
   
+  leng <- length(x)
+  dens <- numeric(leng)
+  
+  for(i in 1:leng){
+  if(x[i]>=0){
   ## On the positive line
   funcp<-function(x,y){y^{delta.minus-1}*(x+y/(gamma.plus+gamma.minus))^{delta.plus-1}*exp(-y)}
   intp<-function(x){integrate(funcp,lower=0,upper=Inf,x=x)$value}
   intvecp<-function(x)sapply(x,intp)
-  densp<-gamma.plus^delta.plus*gamma.minus^delta.minus/((gamma.plus+gamma.minus)^delta.minus*gamma(delta.plus)*gamma(delta.minus))*exp(-gamma.plus*x)*intvecp(x)
-  
+  dens[i]<-gamma.plus^delta.plus*gamma.minus^delta.minus/((gamma.plus+gamma.minus)^delta.minus*gamma(delta.plus)*gamma(delta.minus))*exp(-gamma.plus*x[i])*intvecp(x[i])
+  }else{
   ## On the negative line
   funcm<-function(x,y){y^{delta.plus-1}*(-x+y/(gamma.plus+gamma.minus))^{delta.minus-1}*exp(-y)}
   intm<-function(x){integrate(funcm,lower=0,upper=Inf,x=x)$value}
   intvecm<-function(x)sapply(x,intm)
-  densm<-gamma.plus^delta.plus*gamma.minus^delta.minus/((gamma.plus+gamma.minus)^delta.plus*gamma(delta.minus)*gamma(delta.plus))*exp(gamma.minus*x)*intvecm(x)
-
-  dens<-ifelse(0<=x,densp,densm)
+  dens[i]<-gamma.plus^delta.plus*gamma.minus^delta.minus/((gamma.plus+gamma.minus)^delta.plus*gamma(delta.minus)*gamma(delta.plus))*exp(gamma.minus*x[i])*intvecm(x[i])
+  }
+  }
   dens
-  
 }
 
 ## Generalized inverse Gaussian