Given:

Three positive integers k, m, and n, representing a population containing k+m+n organisms: k individuals are homozygous dominant for a factor, m are heterozygous, and n are homozygous recessive.

Return:

The probability that two randomly selected mating organisms will produce an individual possessing a dominant allele (and thus displaying the dominant phenotype). Assume that any two organisms can mate.

1. sample dataset: k=2,m=2,n=2. Output: 0.78333
2. sample dataset: k=25,m=20,n= 24. Output: 0.7598039215686274

Code: [My code is giving wrong answer for 2. sample dataset]

import numpy as np

GG = np.array([[1,1]]) #so that [00] appears [00] instead of [0,0] matrix
Gg = np.array([[1,0]])
gg = np.array([[0,0]])
each_prob = [] #probability of KK, MM, NN, KM, KN, MN, MK, NK, NM
each_dominant_fraction = []
sum_probs = []

def probability_dom(k,m,n):
    input = [k,m,n]
    for ii in range(len([k,m,n])): #3 represents homo dom, hom recessive, hetero dom for choosing first parent
        for iii in range(len([k,m,n])): #3 represents homo dom, hom recessive, hetero dom for choosing second parent
            t = k+m+n        
            if ii != iii: #if parents are hetero
                probability = (input[ii]/t) * (input[iii]/(t-1))
            else: #if parents are homo
                probability = (input[ii]/t) * ((input[iii]-1)/(t-1))
            each_prob.append(probability)
    print(each_prob)
    print(sum(each_prob)) #should be 1

    for iT in [GG,Gg,gg]:
        for i in [GG, Gg, gg]:
            matrix = iT.T * i #it.T turns it into 2 rows in col and then multiply with another parent, you get offsprings
            Dominant_fraction = np.sum(matrix == 0) / 4
            #np.sum(matrix==0) will count the no of matrix elements that equal 0
            #the whole line of code will define the "Dominant_fraction" variable == to that fraction in the form of a percent
            each_dominant_fraction.append(Dominant_fraction)

    for n in range(0,9):
        sum_probs.append(each_prob[n] * each_dominant_fraction[n])
    print(sum(sum_probs))

probability_dom(25,20,24)