Displaying 2 results from an estimated 2 matches for "integpsv".
2006 Aug 07
1
mathematica -> r (gamma function + integration)
...------------------------------------"];
Print[" "];
Print["p(sv|D_1D_2I) = const. ",N[PSV,6]];
########
# R part
library(fOptions)
###raw values for reproduction
NN <- 58
dd <- 0.411769
lownum <- 20.81512
upnum <- 6.741643
sL <- 0.029
sH <- 0.092
###
integpsv <- function(s) { 1 / (s^NN) * exp(-dd / (2 * s^2)) *
( (igamma((upnum/(2*s^2)),1/2) - igamma(0,1/2) ) +
(igamma((lownum/(2*s^2)),1/2) - igamma(0,1/2) ) )
}
psv <- integrate(integpsv, lower=sL, upper=sH)
PSV <- psv$value / sqrt(2*NN)
print("------------- Results -----------------...
2006 Sep 01
1
integration problem with gamma function
...vbar = NIntegrate[1/(UiA^(Ni/2)) 1/(UiiA^(Nii/2))
Gamma[Ni/2,UiA/(sH^2),UiA/(sL^2)]
Gamma[Nii/2,UiiA/(sH^2),UiiA/(sL^2)],{A,L,H},
MinRecursion->3];
PSVbar = psvbar/(4 Log[sH/sL]);
Print["p(s?v|D_1D_2I) = const. ",N[PSVbar,6]];
</--->
<--- translation to R --->
integpsvbar <- function(A)
{
# Mathematica: gamma[a,z0,z1] = gamma[a,z1] - gamma[a,z0]
UiA <- Ni*(Dsi-2*Di*A+A^2)/2
UiiA <- Nii*(Dsii-2*Dii*A+A^2)/2
1/UiA^(Ni/2) *
1/UiiA^(Nii/2) *
( pgamma(UiA/(sL^2), Ni/2) - pgamma(UiA/(sH^2), Ni/2) ) *
( pgamma(UiiA/(sL^2), Nii/2) - pgamma(UiiA/(sH^2), Nii...