[Rd] rmultinom.c error probability not sum to 1
M.van_Iterson at lumc.nl
M.van_Iterson at lumc.nl
Thu Mar 10 21:25:20 CET 2016
Hi all,
I should have given a better explanation of my problem. Here it is.
I extracted from my code the bit that gives the error. Place this in a file called test.c
#include <math.h>
#include <R.h>
#include <Rmath.h>
#include <float.h>
#include <R_ext/Print.h>
int main(){
double prob[3] = {0.0, 0.0, 0.0};
double prob_tot = 0.;
prob[0] = 0.3*dnorm(2, 0, 1, 0);
prob[1] = 0.5*dnorm(5, 0, 1, 0);
prob[2] = 0.2*dnorm(-3, 0, 1, 0);
//obtain prob_tot
prob_tot = 0.;
for(int j = 0; j < 3; j++)
prob_tot += prob[j];
//normalize probabilities
for(int j = 0; j < 3; j++)
prob[j] = prob[j]/prob_tot;
//or this give the same error
//prob[2] = 1.0 - prob[1] - prob[0];
//checking indeed prob_tot not exactly 1
for(int j = 0; j < 3; j++)
prob_tot += prob[j];
Rprintf("Prob_tot: %f\n", prob_tot);
int rN[3];
rmultinom(1, prob, 1, rN);
return 0;
}
run R CMD SHLIB test.c to generate the test.so. Now from within R
> dyn.load("test.so")
> .C("main")
Prob_tot: 1.017084
Error: rbinom: probability sum should be 1, but is 0.948075
Maybe I miss some trivial C knowledge why this is not exactly one!
Thanks in advance!
Regards,
Maarten
________________________________________
From: peter dalgaard [pdalgd at gmail.com]
Sent: Thursday, March 10, 2016 1:26 PM
To: Iterson, M. van (MOLEPI)
Cc: r-devel at r-project.org
Subject: Re: [Rd] rmultinom.c error probability not sum to 1
On 10 Mar 2016, at 12:47 , M.van_Iterson at lumc.nl wrote:
> Dear all,
>
> I have a questions regarding using the c function rmultinom.c.
>
> I got the following error message "rbinom: probability sum should be 1, but is 0.999264"
>
> Which is thrown by:
>
> if(fabs((double)(p_tot - 1.)) > 1e-7)
> MATHLIB_ERROR(_("rbinom: probability sum should be 1, but is %g"),
> (double) p_tot);
>
> I understand my probabilities do not sum to one close enough. I tried the following,
> p[2] = 1. - p[0] - p[1], where 'p', are the probabilities but this is not sufficient to pass the error message!
>
p[0] ????
> Thanks in advance!
>
> Regards,
> Maarten
>
> (I don't think this is an issue with versions but I used R version 3.2.3 and can provide more details on my linux build if necessary.)
>
>
>
> [[alternative HTML version deleted]]
>
> ______________________________________________
> R-devel at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel
--
Peter Dalgaard, Professor,
Center for Statistics, Copenhagen Business School
Solbjerg Plads 3, 2000 Frederiksberg, Denmark
Phone: (+45)38153501
Office: A 4.23
Email: pd.mes at cbs.dk Priv: PDalgd at gmail.com
More information about the R-devel
mailing list