[prev in list] [next in list] [prev in thread] [next in thread] 

List:       r-devel
Subject:    Re: [Rd] rmultinom.c error probability not sum to 1
From:       <M.van_Iterson () lumc ! nl>
Date:       2016-03-10 21:12:35
Message-ID: 6C970093269C804E91E22510E87295281C3A4647 () MAIL-MB01 ! lumcnet ! prod ! intern
[Download RAW message or body]

Thanks for your replies. Indeed, I had to add the dimension of the probability vector \
like this

 rmultinom(1, prob, 3, rN);

Regards, 
Maarten
________________________________________
From: peter dalgaard [pdalgd@gmail.com]
Sent: Thursday, March 10, 2016 9:45 PM
To: Iterson, M. van (MOLEPI)
Cc: r-devel@r-project.org
Subject: Re: [Rd] rmultinom.c error probability not sum to 1

> On 10 Mar 2016, at 21:25 , M.van_Iterson@lumc.nl wrote:
> 
> 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

Aha. Missing info #1, C not R...

> 
> #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];
> 

prob_tot = 0; missing here

> //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);

Er, where do you tell rmultinom that variates are three-dimensional? It's not going \
to infer it from array sizes.

-pd

> 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@gmail.com]
> Sent: Thursday, March 10, 2016 1:26 PM
> To: Iterson, M. van (MOLEPI)
> Cc: r-devel@r-project.org
> Subject: Re: [Rd] rmultinom.c error probability not sum to 1
> 
> On 10 Mar 2016, at 12:47 , M.van_Iterson@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@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@cbs.dk  Priv: PDalgd@gmail.com
> 

--
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@cbs.dk  Priv: PDalgd@gmail.com










______________________________________________
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel


[prev in list] [next in list] [prev in thread] [next in thread]