We propose a probability distribution for multivariate binary random variables. The probability distribution is expressed as principal minors of the parameter matrix, which is a matrix analogous to the inverse covariance matrix in the multivariate Gaussian distribution. In our model, the partition function, central moments, and the marginal and conditional distributions are expressed analytically. That is, summation over all possible states is not necessary for obtaining the partition function and various expected values, which is a problem with the conventional multivariate Bernoulli distribution. The proposed model has many similarities to the multivariate Gaussian distribution. For example, the marginal and conditional distributions are expressed in terms of the parameter matrix and its inverse matrix, respectively. That is, the inverse matrix represents a sort of partial correlation. The proposed distribution can be derived using Grassmann numbers, anticommuting numbers. Analytical expressions for the marginal and conditional distributions are also useful in generating random numbers for multivariate binary variables. Hence, we investigated sampling distributions of parameter estimates using synthetic datasets. The computational complexity of maximum likelihood estimation from observed data is proportional to the number of unique observed states, not to the number of all possible states as is required in the case of the conventional multivariate Bernoulli distribution. We empirically observed that the sampling distributions of the maximum likelihood estimates appear to be consistent and asymptotically normal.