Buggy python submission (verdict: Time Limit Exceeded):
import sys
input = sys.stdin.readline


def main():
    N = int(input())
    # dp[i][j][k] := 寿司が 1 つある皿が i 皿、2 つが j 皿、
    # 3 つが k 皿の状態から全て食べ終えるまでの操作回数の期待値
    dp = [[[0.0] * (N+1) for _ in range(N+1)] for _ in range(N+1)]

    a = list(map(int, input().split()))
    I, J, K = a.count(1), a.count(2), a.count(3)

    # 小さい方から dp を埋める
    for k in range(N+1):
        for j in range(N+1):
            for i in range(N+1):
                s = i + j + k
                if s != 0 and s <= N:
                    # 式変形はググった…
                    dp[i][j][k] = N
                    if i > 0:
                        dp[i][j][k] += dp[i-1][j][k] * i
                    if j > 0:
                        dp[i][j][k] += dp[i+1][j-1][k] * j
                    if k > 0:
                        dp[i][j][k] += dp[i][j+1][k-1] * k
                    dp[i][j][k] /= s

    print(dp[I][J][K])


if __name__ == '__main__':
    main()


Fixed submission:
import sys
input = sys.stdin.readline

N = int(input())
# dp[i][j][k] := 寿司が 1 つある皿が i 皿、2 つが j 皿、
# 3 つが k 皿の状態から全て食べ終えるまでの操作回数の期待値
dp = [[[float(0)] * (N+1) for _ in range(N+1)] for _ in range(N+1)]

a = list(map(int, input().split()))
I, J, K = a.count(1), a.count(2), a.count(3)

# 小さい方から dp を埋める
for k in range(N+1):
    for j in range(N+1):
        for i in range(N+1):
            if (i+j+k) <= N:
                # 式変形はググった…
                if (i+j+k) != 0:
                    dp[i][j][k] = float(N) / float(i+j+k)
                if i > 0:
                    dp[i][j][k] += dp[i-1][j][k] * float(i) / float(i+j+k)
                if j > 0:
                    dp[i][j][k] += dp[i+1][j-1][k] * float(j) / float(i+j+k)
                if k > 0:
                    dp[i][j][k] += dp[i][j+1][k-1] * float(k) / float(i+j+k)
print(dp[I][J][K])

Buggy python submission (verdict: Time Limit Exceeded):
import sys
from collections import defaultdict
input = sys.stdin.readline


def main():
    N = int(input())
    S = list(map(int, input().split()))
    tmp = [0, 0, 0]  # 皿パターン: 1,2,3個の寿司が置かれている皿の数
    for s in S:
        tmp[s-1] += 1

    d = defaultdict(float)
    d[tuple(tmp)] = 1.0

    e = 0.0  # 求める期待値
    X = sum(S)  # 残りの寿司の数
    while X > 0:
        A = d.items()
        d = defaultdict(float)
        for k, v in A:
            if k == (0, 0, 0):
                continue

            x = v / sum(k)
            e += x * N

            if k[0] > 0:
                d[(k[0]-1, k[1], k[2])] += k[0] * x
            if k[1] > 0:
                d[(k[0]+1, k[1]-1, k[2])] += k[1] * x
            if k[2] > 0:
                d[(k[0], k[1]+1, k[2]-1)] += k[2] * x
        X -= 1

    print(e)


if __name__ == '__main__':
    main()


Fixed submission:
import sys
input = sys.stdin.readline

N = int(input())
# dp[i][j][k] := 寿司が 1 つある皿が i 皿、2 つが j 皿、
# 3 つが k 皿の状態から全て食べ終えるまでの操作回数の期待値
dp = [[[float(0)] * (N+1) for _ in range(N+1)] for _ in range(N+1)]

a = list(map(int, input().split()))
I, J, K = a.count(1), a.count(2), a.count(3)

# 小さい方から dp を埋める
for k in range(N+1):
    for j in range(N+1):
        for i in range(N+1):
            if (i+j+k) <= N:
                # 式変形はググった…
                if (i+j+k) != 0:
                    dp[i][j][k] = float(N) / float(i+j+k)
                if i > 0:
                    dp[i][j][k] += dp[i-1][j][k] * float(i) / float(i+j+k)
                if j > 0:
                    dp[i][j][k] += dp[i+1][j-1][k] * float(j) / float(i+j+k)
                if k > 0:
                    dp[i][j][k] += dp[i][j+1][k-1] * float(k) / float(i+j+k)
print(dp[I][J][K])

Buggy python submission (verdict: Time Limit Exceeded):
import sys
from collections import defaultdict
input = sys.stdin.readline


def main():
    N = int(input())
    S = list(map(int, input().split()))
    tmp = [0, 0, 0]  # 皿パターン: 1,2,3個の寿司が置かれている皿の数
    for s in S:
        tmp[s-1] += 1
    d = defaultdict(int)
    d[tuple(tmp)] = 1

    e = 0  # 求める期待値
    X = sum(S)  # 残りの寿司の数
    while X > 0:
        A = d.items()
        d = defaultdict(int)
        for k, v in A:
            if k == (0, 0, 0):
                continue
            
            e += v / (sum(k)/N)  # 皿パターン数 * (寿司を食べられる期待値 := 寿司がある確率の逆数)

            if k[0] > 0:
                d[(k[0]-1, k[1], k[2])] += k[0] * v / sum(k)
            if k[1] > 0:
                d[(k[0]+1, k[1]-1, k[2])] += k[1] * v / sum(k)
            if k[2] > 0:
                d[(k[0], k[1]+1, k[2]-1)] += k[2] * v / sum(k)
        X -= 1

    print(e)


if __name__ == '__main__':
    main()


Fixed submission:
import sys
input = sys.stdin.readline

N = int(input())
# dp[i][j][k] := 寿司が 1 つある皿が i 皿、2 つが j 皿、
# 3 つが k 皿の状態から全て食べ終えるまでの操作回数の期待値
dp = [[[float(0)] * (N+1) for _ in range(N+1)] for _ in range(N+1)]

a = list(map(int, input().split()))
I, J, K = a.count(1), a.count(2), a.count(3)

# 小さい方から dp を埋める
for k in range(N+1):
    for j in range(N+1):
        for i in range(N+1):
            if (i+j+k) <= N:
                # 式変形はググった…
                if (i+j+k) != 0:
                    dp[i][j][k] = float(N) / float(i+j+k)
                if i > 0:
                    dp[i][j][k] += dp[i-1][j][k] * float(i) / float(i+j+k)
                if j > 0:
                    dp[i][j][k] += dp[i+1][j-1][k] * float(j) / float(i+j+k)
                if k > 0:
                    dp[i][j][k] += dp[i][j+1][k-1] * float(k) / float(i+j+k)
print(dp[I][J][K])