I got confuse when try to calculate Expectation value with Python (with and without replacement)! Sorry if my question is more relate to programing than Math!

A committee of four people is drawn at random from a set of six men and three women. Suppose we are concerned that there may be quite a gender imbalance in the membership of the committee. Toward that end, let $M$ and $W$ denote the numbers of men and women in our committee, and let
the difference be $D = M − W$. Let’s find $E(D)$.

How I calculate:
$D = \{-2,0,2,4\}, D = -2$ means there are $1$ man and $3$ women.
$$E(D) = -2 \times P(D = -2) + 0 + 2 \times P(D = 2) + 4 \times P(D = -4)$$
Without Replacement
$P(D = -2) = \frac{\binom{6}{1}\binom{3}{3}}{\binom{9}{4}}$ and get result $E(D) = \frac{4}{3}$ ~ 1.3333
With Replacement
I use the formula in this link https://en.wikipedia.org/wiki/Multiset#Counting_multisets
$P(D = -2) = \frac{6 \times \binom{3 + 3 -1}{3}}{\binom{9 + 4 -1}{4}}$ and get result $E(D) = \frac{16}{11}\sim1.4545$
But when I Try to check with Python, both will have almost the same results $1,33$.
My method:
I created a list in which has numbers from 1 to 9. Numbers <= 6 refer to Men, otherwise are women.
I will loop n times, each time I generate a sample and calculate $D$ and add to total. Finally, find the mean total / n.
import random 
import numpy as np
pop = []
for i in range(1,10):
    pop.append(i)

n = 300000
total = 0
for _ in range(n):
    a = np.array(np.random.choice(pop, 4)) #with replacemant
    #a = np.array(random.sample(pop, k=4)) #without replacement
    total += sum(a <= 6) - sum(a > 6)
print('Result',total/n)

 A: $$E[D]=E[M-W]=E[M]-E[W]=4\cdot \frac{6}{9}-4\cdot\frac{3}{9}=\frac{4}{3}$$
There was never any need to calculate $\Pr(D=-2)$ or any of the others.  You could have approached purely by linearity of expectation.
$M$ could be thought of as $M_1+M_2+M_3+M_4$ where $M_i$ is the indicator random variable corresponding to if the $i$'th person selected was male or not.  This was how I was able to so quickly say that $E[M]=4\cdot \frac{6}{9}$.  In both the with replacement and without replacement scenarios, $\Pr(M_i=1)=\frac{6}{9}$ seen immediately by noting that each person is just as likely as each other person to have been the $i$'th person selected... the dependence of the events is completely irrelevant and can be ignored when approaching with linearity of expectation.  Both posed problems have an answer of $\frac{4}{3}$
The punchline is going to be in the end... that your $\binom{n+k-1}{k-1}$ arrangements counted with stars and bars are not going to be equally likely scenarios.  That is the source of your error.  The correct value for with replacement should have had a denominator of $9^4$ and for $D=-2$ (one man three women) would have had a numerator of $\binom{4}{1}\cdot 6^1\cdot 3^3$, corresponding to picking which position in the sequence the male was picked, which male it was in that position, and which female(s) they were in those positions.  Note also that for the with replacement problem if you were to insist on approaching the hard and tedious long way, there also exists the possibility of $D=-4$.
A: How does sampling with replacement make any practical
sense when choosing a committee? (I ignore
computations that involve sampling with replacement.)
I agree with
your answer $E(D) = 1.3333.$
Because you're interested in verification via
simulation, I simulated this in R, using 1's
for the men, -1's for the women. (The sample procedure in R samples without replacement by default, treating the nine elements of mem
as distinct; for example, it will never pick
four -1's.)
I get the approximation
$E(D)=1.3322 \pm 0.0030.$ [With a million iterations
one can expect 2 or 3 place accuracy.]
set.seed(1234)
mem = c(1,1,1,1,1,1,-1,-1,-1)
d = replicate(10^6, sum(sample(mem, 4)))
mean(d)
[1] 1.332196     # aprx E(D)
2*sd(d)/1000
[1] 0.002979235  # 95% margin of sim error
table(d)/10^6
d
       -2        0        2        4 
 0.047345 0.358031 0.475805 0.118819 

Notes: (1) Also, your combinatorial value for $P(D=-2) = 0.0476,$
compared with simulated value $0.0473,$ which agrees within the margin of simulation error.
(2) I don't use Python enough to feel
competent commenting on your program.
