https://github.com/cran/nacopula
Tip revision: 5bc804b03d066ef4a2ab9cf3af6f4f2df5bcda4e authored by Martin Maechler on 23 September 2011, 00:00:00 UTC
version 0.7-9-1
version 0.7-9-1
Tip revision: 5bc804b
rSibuya.c
/*
Copyright (C) 2010 Marius Hofert and Martin Maechler
This program is free software; you can redistribute it and/or modify it under
the terms of the GNU General Public License as published by the Free Software
Foundation; either version 3 of the License, or (at your option) any later
version.
This program is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
details.
You should have received a copy of the GNU General Public License along with
this program; if not, see <http://www.gnu.org/licenses/>.
*/
#include <Rmath.h>
#include "nacopula.h"
/**
* Sample V from a Sibuya(alpha) distribution with cdf F(n) = 1-1/(n*B(n,1-alpha)),
* n in IN, with Laplace-Stieltjes transform 1-(1-exp(-t))^alpha via the
* algorithm of Hofert (2011).
* Note: The caller of this function must use GetRNGstate() and PutRNGstate().
* @param alpha parameter theta0/theta1 in (0,1]
* @param gamma_1_a Gamma(1-alpha)
* @return a random variate from F
* @author Marius Hofert, Martin Maechler
*/
double rSibuya(double alpha, double gamma_1_a /**< == Gamma(1 - alpha) */){
/**< FIXME(MM): (for alpha not too close to 1): re-express using 1-U */
double U = unif_rand();
if(U <= alpha)
return 1.;
else { /**< alpha < U < 1 */
const double xMax = 1./DBL_EPSILON; // ==> floor(x) == ceil(x) for x >= xMax
double Ginv = pow((1-U)*gamma_1_a, -1./alpha), fGinv = floor(Ginv);
if(Ginv > xMax) return fGinv; /* else */
if(1-U < 1./(fGinv*beta(fGinv, 1.-alpha))) return ceil(Ginv); /* else */
return fGinv;
}
}
/**
* Sample an n-fold sum of i.i.d. V ~ Sibuya(alpha). Sum-version of rSibuya.
* Note: The caller of this function must use GetRNGstate() and PutRNGstate().
* @param alpha parameter theta0/theta1 in (0,1]
* @param gamma_1_a Gamma(1-alpha)
* @return n-fold sum of i.i.d. V ~ F
* @author Marius Hofert
*/
double rSibuya_sum(int n, double alpha, double gamma_1_a /**< == Gamma(1 - alpha) */){
double n_sum_V = 0.;
for(int i = 0; i < n; i++)
n_sum_V += rSibuya(alpha, gamma_1_a);
return n_sum_V;
}
/**
* Generate a vector of variates from a Sibuya(alpha) distribution.
* @param V vector of random variates from F (result)
* @param n length of the vector V
* @param alpha parameter theta0/theta1 in (0,1]
* @return none
* @author Marius Hofert, Martin Maechler
*/
void rSibuya_vec(double* V, int n, double alpha){
if(n >= 1) {
double gamma_1_a = gammafn(1.-alpha);
GetRNGstate();
for(int i=0; i < n; i++)
V[i] = rSibuya(alpha, gamma_1_a);
PutRNGstate();
}
return;
}
/**
* Generate a vector of variates from a Sibuya(alpha) distribution. Bridge to R.
* @param n sample size
* @param alpha parameter theta0/theta1 in (0,1]
* @return vector of random variates
* @author Martin Maechler
*/
SEXP rSibuya_vec_c(SEXP n_, SEXP alpha_){
int n = asInteger(n_);
double alpha = asReal(alpha_);
SEXP res = PROTECT(allocVector(REALSXP, n));
if(n >= 1) rSibuya_vec(REAL(res), n, alpha);
UNPROTECT(1);
return(res);
}