What's wrong with this simple method to sample from multinomial in C#?

I wanted to implement a simple method to sample from a multinomial distribution in C# (the first argument is an array of integers we want to sample and the second one is the probabilities of selecting each of those integers).

When I do this with numpy in python, the results make sense.

np.random.choice(np.array([1,2,3,4,5,6]),p=np.array([.624,.23,.08,.04, .02, .006]),size=len(b))

I get a lot of 1's (probability 62%), a bunch of 2's, some 3's etc.

However, when I try the implementation below in C# (pretty straightforward inverse transform sampling for multinomial, only relies on a uniform random variable), I get really weird results. For all 1000 samples, I'll often find all 1's. Sometimes, I'll find all 3's (!!??). The results never look like what you would expect (and what you get from the python function - try running it yourself a few times). This is really scary since we rely on these primitives. Does anyone have insight into what might be wrong with the C# version?

    static void Main(string[] args)
    {
        int[] iis = new int[7];
        int[] itms = new int[] { 1, 2, 3, 4, 5, 6 };
        double[] probs = new double[] { .624, .23, .08, .04, .02, .006 };
        for (int i = 0; i < 1000; i++)
        {
            iis[MultinomialSample(itms, probs)] += 1;
        }

        foreach (var ii in iis)
        {
            Console.Write(ii + ",");
        }

        Console.Read();
    }


     private static int MultinomialSample(int[] s, double[] ps)
    {
        double[] cumProbs = new double[ps.Length];
        cumProbs[0] = ps[0];

        for (int i = 1; i < ps.Length; i++)
        {
            cumProbs[i] = cumProbs[i - 1] + ps[i];
        }

        Random random = new Random();
        double u = random.NextDouble();

        for (int i = 0; i < cumProbs.Length - 1; i++)
        {
            if (u < cumProbs[i])
            {
                return s[i];
            }
        }

        return s[s.Length - 1];
    }