It is easy to calculate the probability for making that observation, given the fact the two coins are equal. This can be done by Fishers exact test. Given these observations
$$
\begin{array} {r|c|c}
&\text{coin }1 &\text{coin }2 \\
\hline
\text{heads} &H_1 &H_2\\
\hline
\text{tails} &n_1-H_1 &n_2-H_2\\\end{array}
$$
the probability to observe these numbers while the coins are equal given the number of tries $n_1$, $n_2$ and the total amount of heads $H_1+H_2$ is
$$
p(H_1, H_2|n_1, n_2, H_1+H_2) = \frac{(H_1+H_2)!(n_1+n_2-H_1-H_2)!n_1!n_2!}{H_1!H_2!(n_1-H_1)!(n_2-H_2)!(n_1+n_2)!}.
$$
But what you are asking for is the probability that one coin is better.
Since we argue about a believe on how biased the coins are we have to use a Bayesian approach to calculate the result. Please note, that in Bayesian inference the term belief is modeled as probability and the two terms are used interchangeably (s. Bayesian probability). We call the probability that coin $i$ tosses heads $p_i$. The posterior distribution after observation, for this $p_i$ is given by Bayes' theorem:
$$
f(p_i|H_i,n_i)= \frac{f(H_i|p_i,n_i)f(p_i)}{f(n_i,H_i)}
$$
The probability density function (pdf) $f(H_i|p_i,n_i)$ is given by the Binomial probability, since the individual tries are Bernoulli experiments:
$$
f(H_i|p_i,n_i) = \binom{n_i}{H_i}p_i^{H_i}(1-p_i)^{n_i-H_i}
$$
I assume the prior knowledge on $f(p_i)$ is that $p_i$ could lie anywhere between $0$ and $1$ with equal probability, hence $f(p_i) = 1$. So the nominator is $f(H_i|p_i,n_i)f(p_i)= f(H_i|p_i,n_i)$.
In Order to calculate $f(n_i,H_i)$ we use the fact that the integral over a pdf has to be one $\int_0^1f(p|H_i,n_i)\mathrm dp = 1$. So the denominator will be a constant factor to achieve just that. There is a known pdf that differs from the nominator by only a constant factor, which is the beta distribution. Hence
$$
f(p_i|H_i,n_i) = \frac{1}{B(H_i+1, n_i-H_i+1)}p_i^{H_i}(1-p_i)^{n_i-H_i}.
$$
The pdf for the pair of probabilities of independent coins is
$$
f(p_1,p_2|H_1,n_1,H_2,n_2) = f(p_1|H_1,n_1)f(p_2|H_2,n_2).
$$
Now we need to integrate this over the cases in which $p_1>p_2$ in order to find out how probable coin $1$ is better then coin $2$:
$$\begin{align}
\mathbb P(p_1>p_2)
&= \int_0^1 \int_0^{p‘_1} f(p‘_1,p‘_2|H_1,n_1,H_2,n_2)\mathrm dp‘_2 \mathrm dp‘_1\\
&=\int_0^1 \frac{B(p‘_1;H_2+1,n_2-H_2+1)}{B(H_2+1,n_2-H_2+1)}
f(p‘_1|H_1,n_1)\mathrm dp‘_1
\end{align}$$
I cannot solve this last integral analytically but one can solve it numerically with a computer after plugging in the numbers. $B(\cdot,\cdot)$ is the beta function and $B(\cdot;\cdot,\cdot)$ is the incomplete beta function. Note that $\mathbb P(p_1=p_2) = 0$ because $p_1$ is a continues variable and never exactly the same as $p_2$.
Concerning the prior assumption on $f(p_i)$ and remarks on it: A good alternative to model many believes is to use a beta distribution $Beta(a_i+1,b_i+1)$. This would lead to a final probability
$$
\mathbb P(p_1>p_2)
=\int_0^1 \frac{B(p‘_1;H_2+1+a_2,n_2-H_2+1+b_2)}{B(H_2+1+a_2,n_2-H_2+1+b_2)}
f(p‘_1|H_1+a_1,n_1+a_1+b_1)\mathrm dp‘_1.
$$
That way one could model a strong bias towards regular coins by large but equal $a_i$, $b_i$. It would be equivalent to tossing the coin $a_i+b_i$ additional times and receiving $a_i$ heads hence equivalent to just having more data. $a_i + b_i$ is the amount of tosses we would not have to make if we include this prior.
The OP stated that the two coins are both biased to an unknown degree. So I understood all knowledge has to be inferred from the observations. This is why I opted for an uninformative prior that dose not bias the result e.g. towards regular coins.
All information can be conveyed in the form of $(H_i, n_i)$ per coin. The lack of an informative prior only means more observations are needed to decide which coin is better with high probability.
Here is the code in R that provides a function P(n1, H1, n2, H2)
$=\mathbb P(p_1>p_2)$ using the uniform prior $f(p_i)=1$:
mp <- function(p1, n1, H1, n2, H2) {
f1 <- pbeta(p1, H2 + 1, n2 - H2 + 1)
f2 <- dbeta(p1, H1 + 1, n1 - H1 + 1)
return(f1 * f2)
}
P <- function(n1, H1, n2, H2) {
return(integrate(mp, 0, 1, n1, H1, n2, H2))
}
You can draw $P(p_1>p_2)$ for different experimental results and fixed $n_1$, $n_2$ e.g. $n_1=n_2=4$ with this code sniped:
library(lattice)
n1 <- 4
n2 <- 4
heads <- expand.grid(H1 = 0:n1, H2 = 0:n2)
heads$P <- apply(heads, 1, function(H) P(n1, H[1], n2, H[2])$value)
levelplot(P ~ H1 + H2, heads, main = "P(p1 > p2)")

You may need to install.packages("lattice")
first.
One can see, that even with the uniform prior and a small sample size, the probability or believe that one coin is better can become quite solid, when $H_1$ and $H_2$ differ enough. An even smaller relative difference is needed if $n_1$ and $n_2$ are even larger. Here is a plot for $n_1=100$ and $n_2=200$:

Martijn Weterings suggested to calculate the posterior probability distribution for the difference between $p_1$ and $p_2$. This can be done by integrating the pdf of the pair over the set $S(d)=\{(p_1,p_2)\in[0,1]^2|d=|p_1-p_2|\}$:
$$\begin{align}
f(d|H_1,n_1,H_2,n_2)
&= \int_{S(d)}f(p_1,p_2|H_1,n_1,H_2,n_2) \mathrm d\gamma\\
&= \int_0^{1-d} f(p,p+d|H_1,n_1,H_2,n_2) \mathrm dp + \int_d^1 f(p,p-d|H_1,n_1,H_2,n_2) \mathrm dp\\
\end{align}$$
Again, not an integral I can solve analytically but the R code would be:
d1 <- function(p, d, n1, H1, n2, H2) {
f1 <- dbeta(p, H1 + 1, n1 - H1 + 1)
f2 <- dbeta(p + d, H2 + 1, n2 - H2 + 1)
return(f1 * f2)
}
d2 <- function(p, d, n1, H1, n2, H2) {
f1 <- dbeta(p, H1 + 1, n1 - H1 + 1)
f2 <- dbeta(p - d, H2 + 1, n2 - H2 + 1)
return(f1 * f2)
}
fd <- function(d, n1, H1, n2, H2) {
if (d==1) return(0)
s1 <- integrate(d1, 0, 1-d, d, n1, H1, n2, H2)
s2 <- integrate(d2, d, 1, d, n1, H1, n2, H2)
return(s1$value + s2$value)
}
I plotted $f(d|n_1,H_1,n_2,H_2)$ for $n_1=4$, $H_1=3$, $n_2=4$ and all values of $H_2$:
n1 <- 4
n2 <- 4
H1 <- 3
d <- seq(0, 1, length = 500)
get_f <- function(H2) sapply(d, fd, n1, H1, n2, H2)
dat <- sapply(0:n2, get_f)
matplot(d, dat, type = "l", ylab = "Density",
main = "f(d | 4, 3, 4, H2)")
legend("topright", legend = paste("H2 =", 0:n2),
col = 1:(n2 + 1), pch = "-")

You can calculate the probability of $|p_1-p_2|$ to be above a value $d$ by integrate(fd, d, 1, n1, H1, n2, H2)
. Mind that the double application of the numerical integral comes with some numerical error. E.g. integrate(fd, 0, 1, n1, H1, n2, H2)
should always equal $1$ since $d$ always takes a value between $0$ and $1$. But the result often deviates slightly.