Skip to content

rmultinom.c error probability not sum to 1

6 messages · M.van_Iterson at lumc.nl, Peter Dalgaard, Berend Hasselman

Original
M.van_Iterson at lumc.nl
# 
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!

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.)
Peter Dalgaard
#
On 10 Mar 2016, at 12:47 , M.van_Iterson at lumc.nl wrote:

            

      
      
      
    

        
p[0] ????

            

      
      
      
    

        

      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      M.van_Iterson at lumc.nl
    

    

      
      #
    

  


  

      
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

            

      
      
      
    

        
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

        

      
      
    

  
  


  


      

    

    

      
      

        


  

        

      Berend Hasselman
    

    

      
      #
    

  


  

          

      
      
      
    

        
You haven't initialized prob_tot before starting the for loop.
Set prob_tot to 0.0 before starting the for loop.

Berend

            

      
      
      
    

        

      
      
    

  
  


  


      

    

    

      
      

        


  

        

      Peter Dalgaard
    

    

      
      #
    

  


  

          

      
      
      
    

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

            

      
      
      
    

        
prob_tot = 0; missing here

            

      
      
      
    

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

-pd

            

      
      
      
    

        

      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      M.van_Iterson at lumc.nl
    

    

      
      #
    

  


  

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

 rmultinom(1, prob, 3, rN);

Regards, 
Maarten