/* 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 . */ #include #include "nacopula.h" /** * Sample V ~ F with Laplace-Stieltjes transform * (1-(1-exp(-t)*(1-e^(-theta1)))^alpha)/(1-e^(-theta0)) * via the algorithm of Hofert (2010). C version. * @param p parameter 1-e^(-theta1) * @param alpha parameter theta0/theta1 in (0,1] * @param iAlpha 1-alpha * @param theta0_le_1 in {0,1} with 1 if and only if theta0 <= 1 * @return a random variate from F * @author Marius Hofert, Martin Maechler */ double rFFrank(double p, double alpha, double iAlpha /**< == 1 - alpha */, int theta0_le_1){ double U, V; if(theta0_le_1) { do { U = unif_rand(); V = rLog(p); } while (U*(V-alpha) > 1./beta(V, iAlpha)); } else { double gamma_1_a = gammafn(iAlpha); do { U = unif_rand(); V = rFJoe(alpha, iAlpha, gamma_1_a); } while(U > pow(p, V-1.)); } return V; } /** * Vectorize rFFrank. Generate a vector of variates * V ~ F with Laplace-Stieltjes transform * (1-(1-exp(-t)*(1-e^(-theta1)))^alpha)/(1-e^(-theta0)). * @param V vector of random variates from F (result) * @param n length of the vector V * @param theta_0 parameter theta0 in (0,infinity) * @param theta_1 parameter theta1 in [theta0, infinity) * @return none * @author Marius Hofert, Martin Maechler */ void rFFrank_vec(double *V, const int n, const double theta_0, const double theta_1){ double p = - expm1(-theta_1), alpha = theta_0 / theta_1, iAlpha = (theta_1 - theta_0) / theta_1; int th_0_le_1 = (theta_0 <= theta_1); if(n >= 1) { GetRNGstate(); for(int i=0; i < n; i++) V[i] = rFFrank(p, alpha, iAlpha, th_0_le_1); PutRNGstate(); } return; } /** * Generate a vector of variates V ~ F with Laplace-Stieltjes transform * (1-(1-exp(-t)*(1-e^(-theta1)))^alpha)/(1-e^(-theta0)). * @param n sample size * @param theta_0 parameter theta0 in (0,infinity) * @param theta_1 parameter theta1 in [theta0, infinity) * @return vector of random variates V * @author Martin Maechler */ SEXP rFFrank_c(SEXP n_, SEXP theta_0_, SEXP theta_1_){ int n = asInteger(n_); double theta_0 = asReal(theta_0_), theta_1 = asReal(theta_1_); SEXP res = PROTECT(allocVector(REALSXP, n)); if(n >= 1) rFFrank_vec(REAL(res), n, theta_0, theta_1); UNPROTECT(1); return res; }