Let $A$ and $B$ be real-valued random variables, with $f_A$ and $f_B$ their probability density functions. Let's say we can observe the values of $A$ many times and estimate $f_A$ fairly precisely. We suppose that $P(A\mid B=y)$ follows a binomial distribution $B(N, y)$, where $N$ is a (known or unknown) constant and $y\in [0,1]$ ($B$ does not take values outside $[0,1]$). What I would like to do is to numerically infer $f_B$ from that observation.
From the law of total probability we have $$f_A(x)=\int_0^1 f_A(x\mid B=y)\ f_B(y)\ dy,$$
so I can see this is a kind of inverse problem. I tried discretizing $f_A$ and $f_B$ (approximating them with simple functions with $20$ intervals) using Matlab, but the $20\times 20$ matrix corresponding to $f_A(x\mid B=y)$ often becomes 'ill-conditioned' and Matlab refuses calculation; even when it gives a decent result, $f_B(y)$ has negative values somewhere, which is not desirable.
Are there any good references for this kind of problem? I would appreciate any input on this.