k identical cards from n decks containing d cards I have a physical problem looking quite a basic problem but I can't find a solution.
The problem can be translated as follow:
I have n identical decks containing d cards each.
Taking one random card from each deck, which is the probability to have at least k identical cards?
 A: Hello Daniele and welcome to CV!  While I haven't thought up anything to clean we can quickly get you a some recursion relations that will help you calculate this.  I will consider our draws of ordered hands so the event space has size $d^n$.  We will break this up by the number of card type 1 that appears (for instance the ace of spades).
This card can appear $i$ times in
$$
{n \choose i}(d-1)^{n-i}
$$
ways.  When $i$ is greater than or equal to $k$ these will all satisfy the condition.  Otherwise we will need the other $n-i$ cards to satisfy the condition.  So we can break this into two sums.  The first part is
$$
\sum_{i \geq k} {n \choose i} (d-1)^{n-i}.
$$
The second term is where the recursion comes in.  If we define $f(n,d,k)$ to be the number of ordered hands from $n$ decks of size $d$ with at least $k$ identical cards, the seonc term equals:
$$
\sum_{i<k} {n \choose i}f(n-i, d-1, k).
$$
I'll come back with some python code in a bit.
Edit:  Code Added
from scipy.special import binom
from random import randint
cache = dict()

def f(n, d, k):
    tup = (n, d, k)
    if tup in cache:
        return cache[tup]

    if n < k:
        return 0
    if n == k:
        return d
    if d == 0:
        return 0
    total = 0
    for i in range(n+1):
        if i<k:
            total += binom(n,i) * f(n-i, d-1, k)
        else:
            total += binom(n,i) * (d-1) ** (n-i)

    cache.update({tup: total})
    return total

def g(n, d, k):
    return f(n,d,k) / (d ** n)

def sim(n, d, k, iters=10000):
    y = []
    for _ in range(iters):
        x = []
        for i in range(n):
            x += [randint(0, d-1)]
        y += [most_common(x)]
    z = (np.array(y) >= k)
    return np.mean(z)

A: I have assumed a deck of D cards that are unique.There are N such decks.
Probability of obtaining a specific card is 1/D.
Since there are N decks,the drawing of a card from each deck is considered independent of each other.Each experiment of drawing a card has a probability of 1/D and the experiment is being repeated N times.Hence they are considered Bernoulli trials following Binomial distribution.
From the formula of binomial distribution we have,
$p(X) = ^NC_mP^m(1-P)^{N-m}$
where N is the number of times experiment done,C denotes combination,
m is the number of times of success having probability P,
and X is a successful event.
In our problem P=1/D.
A success event here is getting atleast k identical cards.
Probability of getting k identical cards is given by,
$^NC_k(\frac{1}{D})^k(1-\frac{1}{D})^{N-k}$.
Similarly, 
probability of getting k+1 identical cards is given by replacing k by k+1
and hence we get $^NC_{k+1}(\frac{1}{D})^{k+1}(1-\frac{1}{D})^{N-k-1}$.
We need the probability of at least k identical cards.
This is given by,
P(getting k identical cards)+P(getting k+1 identical cards)+P(getting k+2 identical cards)+...P(getting N identical cards).
m takes values from k to N.
Summation of the 1st equation with range of m from k to N and P=1/D gives the required probability.  
A: I want to put my own answer here, both for my own reference and also for others that stumble upon this question. I have been thinking about this for a while now... and it was fun! The reason why this is an important question and why it is of interest to look into this is because this problem is tightly related to studying the properties of resampling techniques, such as bootstrapping.
You can think about your problem in a way of buckets and balls. I.e. you have $d$ buckets (which represent the different cards) and you can throw $n$ balls (balls represent the draw of a card) into the buckets.
The number of balls that fall into bucket $i$ is a random variable $X_i$. Now you have a random vector $\mathbf{Z}=(X_1,X_2,...,X_d)$ which has the property that $\sum_{i=1}^d X_i = n$. This random vector follows a uniform multinomial distribution. Uniform means that there is equal probability of throwing a ball into each bucket, i.e. $\frac{1}{d}$. Now we can define another random variable $M=\max(\mathbf{Z})$.

Your question can now be reformulated as, how can we find the probability mass function of $M$?

After laying under fur for a couple of days I have not yet managed to find an explicit formula for the probability mass function, and I was starting to think that it may not even exist.
This MO question asks exactly about this reformulation of the question which refers to this paper, which is essentially an algorithm to calculate this, but not an explicit formula.
Even calculating the expectation of $M$ is only achieved via bounds. So there do exist bounds that might give you rough estimates but getting an explicit formula is a significantly harder problem.
The solution presented by @jlimahaverford is an algorithm which is created via a recurrence relation. It is not any simple recurrence relation, so I am not sure if that can be reduced to an explicit formula, but that is probably a good start if someone wants to try it. After my search, I think that his solution is one of the more elegant ways I have seen to do this. So hats of to him!
But! There is a way to get the exact distribution and it is presented in this paper.
