I am new in bayesian analysis and I want to use Approximate Bayesian computation in order to convert an odd giving to me by a bookmaker to a probability that the event occurs. Here is the Python code I created to do it but my model doesn't converge. I guess something is wrong somewhere.
To begin, I have a dataset like this one:
odd outcome 1.47 1 3.21 0 10.22 0 1.17 1 1.17 0 2.12 1 1.82 0
I assume my model is a function with two parameters $a$ and $b$: $proba = f(odd) = b + a/odd$. So, the objective of my ABC is to find the distribution of $a$ and $b$. My prior about these parameters is $a$ and $b$ are independents and $a$ is following a normal distribution with a mean of 1.0 and a standard deviation of 0.1 ($a \sim N(1.0,0.1)$) and b is also following a normal distribution with a mean of 0.0 and a standard deviation of 0.1 ($b \sim N(0.0,0.1)$).
Before building my model with the real dataset, I chose to build one on my own to be sure my model is correct. So, I created my $f$ function with the fixed parameters $a=0.97$ and $b=-0.02$. Then, I generated my dataset: I picked the inverse of a number between 0.01 and 0.99 from an uniform random variable to generate my odd using my $f$ function and I simulated a Bernoulli trial with parameter $p = f(odd)$.
Now, it's time to implement ABC. The idea is to generate $n$ tuples $(a,b)$ from my prior distribution, feed the $f$ function with a subset of $k$ rows from my dataset, perform a Bernoulli trial with $p=f(odd)$ and reject all tuples who have a different result than my dataset.
Here is an example with $n=100000$ and $k=5$:
- I generate a new dataset of 5 rows following the steps described above:
odd outcome 1.47 1 3.21 0 10.22 0 1.17 1 1.17 0
- I generate 100000 tuples $(a,b)$ following my prior distribution:
a b 0.93 0.09 0.86 0.04 1.21 -0.17 0.90 0.07 0.95 -0.16
- For each tuple, I use the odds from my dataset and the value of my parameters $a$ and $b$ to generate a probability using my $f$ function. Then, I perform a Bernoulli trial to get my outcome. For instance, for the tuple (0.93,0.09), I get:
odd f(odd) outcome 1.47 0.72 1 3.21 0.38 0 10.22 0.18 0 1.17 0.88 0 1.17 0.88 1
I should add if $f(odd)$ is greater than 1, I force the result to 1. Same way if the result is lower than 0.
Because the outcome is different than my dataset, I reject that tuple.
I repeat step 3 for all tuples and I get my posterior corresponding to all my tuples I didn't reject.
Because I want to perform these steps $m$ times ($m=1000$ in my code), I have to be able to generate a new list of tuples from my posterior. To do so, I decided to make a kernel-density estimate with a Gaussian kernel to create a new distribution from my posterior.
By this way, I hope the mean of my posterior will converge to $(0.9,-0.02)$ (the value I fixed for $a$ and $b$) but, as I said, it doesn't converge. Here is the evolution of the mean of $a$ and $b$ for each step of my algorithm:
The shaded area represents the mean +/- 1.96 standard deviation. We can also see that the standard deviation is not reducing.
Is something wrong with my methodology?
Edit: I tried to change my prior to one with an incoherent mean and a huge standard deviation and I got these results: We can see that a coherent value for $a$ and $b$ is found really quickly but the algorithm still can't converge to a precise value for these parameters. I guess my algorithm is correct but my model (my $f$ function) can't be more precise. What do you think?