/* 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 a Log(p) distribution with the algorithm "LK" of Kemp (1981). * Note: The caller of this function must use GetRNGstate() and PutRNGstate(). * @param p in (0,1) * @param Ip = 1 - p_ (possibly more accurate) * @return a random variate from Log(p) * @author Marius Hofert, Martin Maechler */ double rLog(double p, double Ip) { if(p <= 0. || p > 1.) { error("rLog(): p must be inside (0,1)"); return -1.; /**< -Wall */ } else if(Ip <= 0. || Ip >= 1.) { error("rLog(): Ip must be inside (0,1)"); return -1.; /**< -Wall */ } else { double U=unif_rand(); if(U > p) { return 1.; } else { double Q, logQ; if(p <= 0.5) { Q = - expm1(log1p(- p) * unif_rand()); /* = 1-(1-p)^unif */ /** * == 1. - exp(log1p(- p) * unif_rand()) * == 1. - pow(1. - p, unif_rand()) */ logQ = log(Q); } else { // p > 0.5 <==> Ip < 0.5 double iQ = pow(Ip, unif_rand()); /* = (1-p)^unif */ Q = 1. - iQ; logQ = log1p(-iQ); } return(U < Q*Q ? floor(1. + log(U)/logQ) : ((U > Q) ? 1. : 2.)); } } } /** * Generate a vector of variates from a Log(p) distribution with the algorithm * "LK" of Kemp (1981). * @param n_ sample size * @param p_ parameter p in (0,1) * @param Ip_ = 1 - p_ (possibly more accurate) * @return vector of random variates from Log(p) * @author Martin Maechler */ SEXP rLog_vec_c(SEXP n_, SEXP p_, SEXP Ip_) { int n = asInteger(n_); double p = asReal(p_),Ip = asReal(Ip_); SEXP res = PROTECT(allocVector(REALSXP, n)); double* X = REAL(res); GetRNGstate(); for(int i=0; i < n; i++) X[i] = rLog(p, Ip); PutRNGstate(); UNPROTECT(1); return res; }