In this game, four six-sided die is being rolled.
What is the probability of where all the dice can be paired off so that one pair adds up to four and the other pair adds up to nine?

First, we have to import numpy in to allow the simulation to work

import numpy as np

Next, we have to create four dice as D1, D2, D3, D4. And with each die having numbers 1-6

D1= np.random.randint(1,7)
D2=np.random.randint(1,7)
D3=np.random.randint(1,7)
D4=np.random.randint(1,7)

After that, we can make the code that runs the simulation multiple time and when there is a run that has one of the wanted result, it return a 1, likewise, it returns 0.

def dice_rolls():
    D1= np.random.randint(1,7)
    D2=np.random.randint(1,7)
    D3=np.random.randint(1,7)
    D4=np.random.randint(1,7)
    
    if (D1+D2==7 and D3+D4==9) or (D1+D2 ==9 and D3+D4==7):
        return 1

    elif (D2+D3==7 and D1+D4==9) or (D2+D3 ==9 and D1+D4==7):
        return 1

    elif (D2+D4==7 and D1+D3==9) or (D2+D4 ==9 and D1+D3==7):
        return 1

    else:
        return 0

Lastly, we can run the simulation multiple times, and add up all the 1 and divide it by the total run to obtain the probabilty. Which in this case we run it 10000 times

sum([dice_rolls() for i in range(100000)])/100000

Based on our simulation, in which we rolled it 100000 times,  if we were to roll four dice, the probabilty of getting a result where we can one pair up the dice with one pair adds up to seven and another pair add up to 9 is 7.399%

We know for 4 dies , total faces could estimate will be 6*6*6*6 . but as the requests, if we pick out D1,D2 , what the face that satisfy the sum(D1,D2)==7 ?

for D1 in range(1,7):
    for D2 in range(1,7):
        if D1+D2==9:
            print((D1,D2))

then another D3 and D4 , same should satisfy the sum(D3,D4)=9.

for D3 in range(1,7):
    for D4 in range(1,7):
        if D3+D4==7:
            print((D3,D4))

total 4 dice (D1,D2,D3,D4), so the sum(faces of dice upon the requests) satisfy the condition that sum(2 dies) is 7 and another sum(2 dice)=9  will be 4(means 4 dies)*6(sum(2 dies)=7)*4(sum(another 2 dies))=96 
sum(faces of 4 dies)=6*6*6*6
thus , the probability=sum(faces of requests dies)/sum(faces of 4 dies )

6*4*4/6**4

Therefore, theoretically speaking, then probability of getting one pair add up to 7 and the other pair add up to 9 is 7.407%, using 4 6-sided dice.

Comparing the the experimental python simulation result to the theoretical result, the two probabilities are relatively close. Since, the probabilty in python is calculated through simulations, the result will change each time. And the more run you do in python, the closer the probability should reach 7.407%.