[Rd] Calling C implementations of rnorm and friends
Luis Usier
luis.henrique.usier at gmail.com
Sat Jul 2 02:39:49 CEST 2016
Gabriel,
That's exactly what I have been doing, and it works fine. However, I just
wanted to *understand* how the random numbers are generated, for no other
reason than to satiate my curiosity.
The one thing that's not very elegant about this way of doing things is
that I have to hard-code special rules for parameter values that are
infinite. I guess there's no other way though; I finally found something
closer to the source code for rnorm which for some reason is *not* in the
stats package folder but in "nmath/rnorm.c":
double rnorm(double mu, double sigma)
{
if (ISNAN(mu) || !R_FINITE(sigma) || sigma < 0.)
ML_ERR_return_NAN;
if (sigma == 0. || !R_FINITE(mu))
return mu; /* includes mu = +/- Inf with finite sigma */
else
return mu + sigma * norm_rand();
}
However, this just scales the result from calling the function norm_rand(),
which I cannot find anywhere : (
Following Joshua's suggestion should only show me functions from the base
package, that is, functions called with .Primitive() and .Internal(). Since
I'm looking at a function from the stats package called with .Call() I
thought it shouldn't be names.c. BUT, it turns out that it *is* actually
there and is also a function that can be called with .Internal() from the
base package. So, surprisingly enough, this:
.Internal(rnorm(1, 0, 1))
actually works, even in the command line, as all functions from base are
always accessible : ) Am I allowed to do this? There's a note in the source
code that says "Can remove all the [dpqr]xxx once the compiler knows how to
optimize
to .External.", whatever that means; I feel like R CMD check/CRAN may not
like it if I call norm and friends from base.
Anyway, thank you both for your help! Any extra nuggets of wisdom would be
appreciated, but this was already very helpful.
Luis
On Sat, Jul 2, 2016 at 12:24 AM, Gabriel Becker <gmbecker at ucdavis.edu>
wrote:
> Well,
>
> For this particular use case why not just transform the parameters at the
> R level and then call the existing function? Is there not a closed form
> mapping?
>
> ~G
> On Jul 1, 2016 2:50 PM, "Joshua Ulrich" <josh.m.ulrich at gmail.com> wrote:
>
>> On Fri, Jul 1, 2016 at 6:13 AM, Luis Usier
>> <luis.henrique.usier at gmail.com> wrote:
>> > Gabriel,
>> >
>> > Thanks for that! I guess I really should have figured that one out
>> sooner,
>> > huh?
>> >
>> > I understand why that wouldn't be CRAN-compliant. But then, what *is*
>> the
>> > proper way to do it? Is there any way I can call unexported functions
>> from
>> > another package and have it accepted by CRAN?
>> >
>> There may be ways to call unexported functions that R CMD check will
>> not detect, but you should not look for ways to violate the spirit of
>> the policy. Work-arounds include asking the package maintainer to
>> export the functionality you want, or you can copy the desired
>> functionality to your package (with attribution, of course).
>>
>> In this particular case, you can probably get what you want via the
>> Rmath library. It provides C-level access to all the distribution
>> functions (and more). You would need something like this in your C
>> code:
>>
>> #include <Rmath.h>
>> GetRNGstate();
>> double rn = rnorm(0.0, 1.0);
>> PutRNGstate()
>>
>> Once Dirk recovers from useR!, he'll tell you how this is all
>> super-easy if you just use Rcpp. ;)
>>
>> > Also, if I instead re-write the random variable generating functions, do
>> > you have any idea of where the source code is in the stats package? As I
>> > said above, I can't seem to find the source code for the functional
>> forms.
>> >
>> See the sections "Functions that call compiled code" and "Compiled
>> code in a base package" from this answer:
>> http://stackoverflow.com/a/19226817/271616
>>
>> That should give you a few pointers on where/how to look.
>>
>> > Thanks,
>> >
>> > Luis
>> >
>> > On Thu, Jun 30, 2016 at 10:38 PM, Gabriel Becker <gmbecker at ucdavis.edu>
>> > wrote:
>> >
>> >> Luis,
>> >>
>> >> C_rnorm is a symbol but it's not exported. This means that you *can*
>> do
>> >> this by using stats:::C_rnorm.
>> >>
>> >> That said, it's not exported, which means that it's not supported to do
>> >> this. So your package likely would not be allowed on CRAN, for example.
>> >>
>> >> Best,
>> >> ~G
>> >> On Jun 30, 2016 2:08 PM, "Luis Usier" <luis.henrique.usier at gmail.com>
>> >> wrote:
>> >>
>> >>> Hi all,
>> >>>
>> >>> Looking at the body for the function rnorm, I see that the body of the
>> >>> function is:
>> >>>
>> >>> .Call(C_rnorm, n, mean, sd)
>> >>>
>> >>> I want to implement functions that generate normal (and other) random
>> >>> variables. Now, I understand that I can perfectly well just call the R
>> >>> wrapper for these functions and that will be almost indistinguishable
>> for
>> >>> most purposes, but for whatever reason I wanted to try and call the C
>> >>> function directly. My first instinct was to call them as:
>> >>>
>> >>> .Call(C_rnorm, 1, 1, 1)
>> >>>
>> >>> This doesn't work because R assumes C_rnorm is an object. Looking at
>> the
>> >>> documentation for .Call, I try passing it in as a string:
>> >>>
>> >>> .Call("C_rnorm", 1, 1, 1, PACKAGE = "stats")
>> >>>
>> >>> This doesn't work either. The help page links to
>> getNativeSymbolInfo(),
>> >>> which I can't make work either. It also refers me to the dyn.load()
>> >>> function and the "Writing R Extensions" manual.
>> >>>
>> >>> After reading and trying to digest those, I try
>> >>>
>> >>> getDLLRegisteredRoutines("stats")
>> >>>
>> >>> which shows me all the C and Fortran functions as registered routines,
>> >>> along with their number of parameters. I retrieve rnorm from the list
>> and
>> >>> pass it on to .Call, which then works fine.
>> >>>
>> >>> However, is there an easier way to do this? Specifically, I would
>> like to
>> >>> call the DLL registered routines from within functions in a package I
>> am
>> >>> writing. The manual tells me I should use useDynLib(). So if I added
>> >>> useDynLib("stats") to my namespace, would that work? Could I then
>> have a
>> >>> function such as:
>> >>>
>> >>> function(x, y, z) .Call(C_rnorm, x, y, z)
>> >>>
>> >>> in my package? If not, what is the proper way of calling these
>> functions
>> >>> from other packages? Should I use "C_rnorm" or "norm"?
>> >>>
>> >>> Also, I was looking for the C source code of rnorm, because I wanted
>> to
>> >>> understand how the function works. Looking at Winston Chang's github R
>> >>> mirror, I found rnorm in the random.c file in the stats package.
>> However,
>> >>> the code I find for it:
>> >>>
>> >>>
>> >>>
>> >>> #define DEFRAND2_REAL(name) \
>> >>> SEXP do_##name(SEXP sn, SEXP sa, SEXP sb) { \
>> >>> return random2(sn, sa, sb, name, REALSXP); \
>> >>> }
>> >>> DEFRAND2_REAL(rnorm)
>> >>>
>> >>>
>> >>> Doesn't help me at all in understanding how it works. It should
>> create a
>> >>> function random2(sn, sa, sb, norm, REALSXP); I understand that is a
>> >>> version
>> >>> of the random2 function that returns a real S expression taking sn,
>> sa and
>> >>> sb as parameters. But how does find the actual functional form for the
>> >>> normal distribution?
>> >>>
>> >>> I am asking because I would like to rewrite some of the other
>> functions,
>> >>> such as parameterizing rbeta by the mean and sample size rather than
>> by
>> >>> the
>> >>> number of successes and failures and rgamma by the mean and total time
>> >>> elapsed instead of the number of events. Once I understand how the C
>> >>> source
>> >>> code works, it would be hopefully not very difficult to reparameterize
>> >>> them.
>> >>>
>> >>> Thanks,
>> >>>
>> >>> Luis Usier
>> >>>
>> >>> [[alternative HTML version deleted]]
>> >>>
>> >>> ______________________________________________
>> >>> R-devel at r-project.org mailing list
>> >>> https://stat.ethz.ch/mailman/listinfo/r-devel
>> >>>
>> >>
>> >
>> > [[alternative HTML version deleted]]
>> >
>> > ______________________________________________
>> > R-devel at r-project.org mailing list
>> > https://stat.ethz.ch/mailman/listinfo/r-devel
>>
>>
>>
>> --
>> Joshua Ulrich | about.me/joshuaulrich
>> FOSS Trading | www.fosstrading.com
>> R/Finance 2016 | www.rinfinance.com
>>
>
[[alternative HTML version deleted]]
More information about the R-devel
mailing list