Solved – rbinom() produces NA values. What’s wrong

binomial distributionmissing datarrandom-generation

I have to generate a random variable from a binomial distribution with R. I have a vector for n and a vector for p and for each component I have to generate a random variable. My idea was the following for the code was the following:

X <- vector(mode="numeric")
for (i in 1:25)
{  
    X[i] <- rbinom(1, n[i], p[i])
}

It seems that my thinking is correct but the vector produced, X, contains only NA. I don't know what the problem is. Maybe I have a problem in the dataset? The vectors n and p are the following:

n <- c(0.8888889, 1.6363636, 0.3333333, 0.6153846, 0.4000000, 0.0625000, 0.2857143,
       0.6153846, 0.0000000, 0.2857143, 0.3333333, 0.1250000, 0.8888889, 0.3333333, 
       0.3333333, 0.0000000, 0.0000000, 0.0000000, 0.6000000, 0.1250000, 0.4000000, 
       0.4000000, 0.0000000, 0.6000000, 1.0000000) 
p <- c(0.8019377, 0.7471281, 0.7673853, 0.7302029, 0.6915447, 0.8211286, 0.7725188, 
       0.7077241, 0.7957015, 0.7893677, 0.7585554, 0.7749553, 0.7570536, 0.6159425, 
       0.7717002, 0.7462224, 0.7349991, 0.7486309, 0.7357340, 0.6850625, 0.7468344, 
       0.7967847, 0.6738766, 0.7122724, 0.7791634)

Best Answer

R expects the second argument of rbinom, size, to be an integer, in accordance with the definition of the binomial distribution. So using a number like 0.9 for size produces NA.

Incidentally, your first code block can be written in one line as X <- rbinom(25, n, p).

Related Solutions

Solved – Generating sorted pseudo-random numbers in Stata

The help for set_seed states

The sequences these functions produce are determined by the seed, which is just a number and which is set to 123456789 every time Stata is launched.

Stata's philosophy emphasizes reproducibility, so this consistency is not surprising. Of course you can set the seed yourself. See the help page for more information.

One way to sort a column separately from all others is to preserve your data, keep only the column to sort, sort it, save the results in a temporary file, restore your data, and merge the temporary file:

gen y = rnormal()
preserve
keep y
sort y
tempfile out
save `out'
restore
merge 1:1 _n using `out', nogen

Solved – Generating Vectors under Constraints on 1 and 2 norm

OK, I think I know how to do it. I am going to stick to your original question but the revised question can be answered using the same method. I'll come back to it at the end.

What is desired is some numbers $x_1, x_2, \ldots x_n$ such that $$x_1 + x_2 + \ldots + x_n = \lambda$$ and $$x_1^2 + x_2^2 + \ldots x_n^2 = \phi$$ and $x_i >0$ for all $i$. I assume you also want to draw uniformly from the set of such points $(x_1, x_2, \ldots, x_n)$ in the sense that every such point has an equal chance of being chosen. I am going to stick to the case $n=3$ for now, because the geometric intuition is easier, but I think the same argument works for $n>3$.

In the $n=3$ case, $\sum x_i^2 = \phi$ defines a sphere of radius $\sqrt{\phi}$ and $\sum x_i = \lambda$ defines a plane. The plane will cut the sphere in a circle, like the circles in this picture (from themathlab.com):

enter image description here

In general, it will look like the "small circle". You want to choose a point uniformly at random from this circle. In the $n$--dimensional case, the equation $\sum x_i = \lambda$ defines a hyperplane and the intersection of this hyperplane with the (hyper)sphere is a (hyper)sphere of one lower dimension. Once you have a description of the hypersphere of intersection, you could use any other method of producing a random point on a sphere.

You can do it by calculating using $n$--dimensional geometry, but the quickest way is probably the following.

Generate a random point from the surface of the unit sphere in dimension $n-1$. This question explains how to do it. Call this point $z$.
Add an extra coordinate to form the point $(0,z)$ in $n$--dimensional space.
Let $U$ be any $n \times n$ orthogonal matrix whose first row consists of $1$s. You could take the real part of a DFT matrix.
Let $$r = \sqrt{\phi -\lambda^2/n}$$ Note that the problem does not have a solution if $\phi -\lambda^2/n$ is negative. Let $v = rU^{-1}(0,z)^T$. Then $||v||^2 = r^2$ and the sum of the entries of $v$ is zero, because it is the first entry of the vector $Uv$, which is zero by construction.
Finally, let $w = v +\frac{\lambda}{n}(1,1, \ldots, 1)$. Then $w$ is a uniform random draw from the required hypersphere. You can check that $\sum w_i^2=\phi$ and $\sum w_i =\lambda$.

You can check that this is satisfactory in the case $n=3$ by using R.

phi <- 0.28
lambda <- 0.9

U <- matrix(c(1,1,1,2,-1,-1,0,-1,1), byrow=T, ncol=3)
U <- U/sqrt(rowSums(U^2))

generate.point <- function(lambda, phi, n){
    z <- rnorm(n-1)
    z <- z/sqrt(sum(z^2))
    z <- c(0,z)
    r <- sqrt(phi - lambda^2/n)
    v <- solve(U, z)*r
    v + (lambda/n)*rep(1,n)
    }

X <- matrix(0, 100, 3)
for (i in 1:100) X[i,] <- generate.point(lambda, phi, 3)
require(lattice)
cloud(X[,3] ~ X[,1] + X[,2])

enter image description here

If you repeat it several times, you should see nice random-looking points on the circle.

As for the version with $\alpha$ and $\beta$, I am not sure whether the constraint $1 > x_i > \alpha$ will be satisfied automatically with $\lambda=1-\alpha$ and $\phi = (\beta-\alpha)\alpha$. If it's not, then you can generate points one at a time and throw away the ones which do not satisfy the extra condition.

Best Answer

Related Solutions

Solved – Generating sorted pseudo-random numbers in Stata

Solved – Generating Vectors under Constraints on 1 and 2 norm

Related Question