2

I want to simulate the pdf of the random variable given in the answer here: Distribution of $n\choose x$ . How can I simulation the pdf of $\binom{n}{X}$ in Matlab in way to confirm that it's indeed what the answer says?

Community
  • 1
noob
  • 639
  • 3
    The sum of two pdfs cannot be a pdf, since it would have total probability 2 (neither can you sum two distributions and get a distribution out). Do you instead mean that you want the pdf of a random variable that's the sum of a Poisson random variable and an exponential random variable (which would suggest you're after the convolution - not the sum - of their pdfs?) --- please clarify your question – Glen_b Mar 06 '17 at 11:37
  • 2
    [At the same time, your Poisson + Exponential situation sounds fascinating; if you're able to give more details about how this situation arises, I expect that could be quite interesting background.] – Glen_b Mar 06 '17 at 11:48
  • @Glen_b, The exponential and poisson are just an example, I will edit my post with what I got – noob Mar 06 '17 at 11:56

2 Answers2

3

You have a binomial random variable $X \sim B(n, p)$, where $n$ and $p$ are given parameters. Now, you want to consider another random variable:

$$Y = \binom{n}{X}$$

$Y$ is a function of $X$ so, to sample from $Y$, first draw samples of $X$. Sampling from a binomial distribution is straightforward. Since you're using Matlab, you can use binornd() with your chosen parameters. Then, pass these samples through the nchoosek() function to obtain samples of $Y$. That is, given samples $\{x_1, \dots, x_n\}$, corresponding samples $\{y_1, \dots, y_n\}$ are given by:

$$y_i = \frac{n!}{x_i! (n-x_i)!}$$

Once you have a large number of samples of $Y$, you can estimate its distribution using the frequency of each value, then compare this to your analytical solution. You can also run goodness-of-fit tests to determine whether the samples match your analytical solution.

Note that $\binom{n}{k}$ grows very quickly (past the computer's capacity to represent numbers using standard formats), so you won't be able to compute it for very large values without playing some special tricks. Otherwise you can stick to moderate values of $n$.

user20160
  • 32,439
  • 3
  • 76
  • 112
  • Could you please explain how do I estimate the pdf using the frequency of each value? – noob Mar 06 '17 at 14:42
  • 1
    Y is discrete, so you can estimate its pmf by counting the number of samples that have each integer value, then dividing by the total number of samples. – user20160 Mar 06 '17 at 14:45
  • Isn't there any way to do this using a function or command, cause I took the number of trials to be 10 and I simulted a sample of length 100. – noob Mar 06 '17 at 14:51
  • 1
    Yes, you can use accumarray() to do the counting. I'd recommend a much higher number of samples (at least 10k, but it's so computationally cheap to do this, you could easily use 100k or more). – user20160 Mar 06 '17 at 15:34
  • I'm sorry but I got confused. I used the function pmf posted here: http://stackoverflow.com/questions/14703240/from-data-to-probability-mass-function-matlab/14703303#14703303 and I got the frequencies, now what? t(he function outputs a vector of the frequency of the unique values not of the data). – noob Mar 06 '17 at 18:24
  • If I understood correctly, your question was about how to approximate the pmf of $Y = \binom{n}{X}$ by numerical simulation, so you can compare it to some pmf that was derived analytically. The pmf gives the probability associated with each of the possible, discrete values that $Y$ can take. This is what you just computed. The unique values are the possible values $Y$ can take, and the frequencies are their probabilities. All other values have probability 0. Now you can compare this to your analytical pmf. – user20160 Mar 06 '17 at 18:51
  • Let us continue this discussion in chat. – noob Mar 06 '17 at 19:36
-1

If you toss a coin $n$ times, and it has a $0.5$ probability head(and tail), then the probability that you get $x$ heads is proportional to $\binom{n}{X}$.

This means that you need to call rand() to generate $n$ uniform random variables, and count how many of them are larger than $0.5$ . This would correspond to a random variable with distribution $P(X=x) = {\binom{n}{x}} / 2^n$.

Mortezaaa
  • 259
  • I am afraid you confuse the pdf equal to ${n \choose x}$ with the pdf of the variable ${n \choose X}$. – Xi'an Mar 09 '17 at 14:33