R devel - Mar 2016 - rmultinom.c error probability not sum to 1

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

> On 10 Mar 2016, at 21:25, M.van_Iterson at 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
> 
> #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];
> 
You haven't initialized prob_tot before starting the for loop.
Set prob_tot to 0.0 before starting the for loop.

Berend
>  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
> 
> 
> ______________________________________________
> R-devel at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel

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

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 at gmail.com]
Sent: Thursday, March 10, 2016 9:45 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 21:25 , M.van_Iterson at 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 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
>
--
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