Initial commit #52
Initial commit #52
Conversation
|
Thanks, looks great! I will make sure to read it and give feedback! |
|
|
||
| Gamma = object | ||
| k: float | ||
| sigma: float |
There was a problem hiding this comment.
Should sigma be theta for consistency with the equations?
| Now that we have prior probability density functions for each model parameter | ||
| that we are estimating, we need to be able to compute the joint prior probability | ||
| for all of the parameters of the model. | ||
| We will actually be using the $ln$ of the probabilities to |
There was a problem hiding this comment.
| We will actually be using the $ln$ of the probabilities to | |
| We will actually be using the $\ln$ of the probabilities to |
| 3. Propose a new value for one of the model parameters by randomly drawing from | ||
| a symetrical distribition centered on the present value. We will use a normal | ||
| distribution. | ||
| 4. Compute the unnormalized posterior probabity density for the proposed model |
There was a problem hiding this comment.
For some reason, this list isn't rendered properly 🤔
There was a problem hiding this comment.
I fear it might be the trailing whitespaces (which oddly enough are meaningful in markdown)
| ggsave(fmt"{prefix}sd.png") | ||
|
|
||
| plotTraces(@[samples, samples2], "images/trace-") | ||
|
|
There was a problem hiding this comment.
I would add a nbClearOutput() here to suppress the ggplotnim echos.
| nbClearOutput() | |
| ggsave(fmt"{prefix}sd.png") | ||
|
|
||
| plotHistograms(@[burnin1, burnin2], "images/hist-") | ||
|
|
There was a problem hiding this comment.
| nbClearOutput() | |
|
|
||
| nbCode: | ||
| plotTraces(@[st_samples_trans1, st_samples_trans2], "images/trace-st-") | ||
|
|
There was a problem hiding this comment.
| nbClearOutput() | |
|
|
||
| nbCode: | ||
| plotHistograms(@[st_burnin1, st_burnin2], "images/hist-st-") | ||
|
|
There was a problem hiding this comment.
| nbClearOutput() | |
| Then to compute the likelihood for a given set of $\beta_{0}$, $\beta_{1}$, $\tau$ | ||
| parameters and data values $x_{i}$ and $y_{i}$ we use the normal probability | ||
| density function which we wrote before to compute our prior probabilities. | ||
| We will work with the $ln$ of the likelihood as we did with the priors. |
There was a problem hiding this comment.
| We will work with the $ln$ of the likelihood as we did with the priors. | |
| We will work with the $\ln$ of the likelihood as we did with the priors. |
|
Oh, I forgot to review this, sorry for the delay, will give a quick look now. |
There was a problem hiding this comment.
Hi, thanks a lot for writing this! I was always interested to learn more about bayesian inference and this definitely is a good starting point. Also it is very nice to see done from scratch in nim stuff that usually is hidden in libraries. I added a few minor points on top of those found by Hugo and marked a few places of confusion for me. In general it could benefit for a bit more external references so the curious reader can better understand the material, which is rather difficult. Apart from that I think it is very well written!
| nbText: hlMd""" | ||
| # Bayesian Inference with Linear Model | ||
|
|
||
| At the time of this writing, Nim does not have any libraries for |
There was a problem hiding this comment.
| At the time of this writing, Nim does not have any libraries for | |
| At the time of this writing (February 2023), Nim does not have any libraries for |
or something similar? let's just it gets out of date soon :)
| the different parts needed to perform Bayesian inference with a linear model. | ||
| We will assume that you have some basic understanding of Bayesian | ||
| inference already. There are many excellent introductions available in books | ||
| and online. |
There was a problem hiding this comment.
It would be nice if you can add a few selected references. At the very minimum I would add a link to https://en.wikipedia.org/wiki/Bayesian_inference for first occurence of Bayesian inference.
Also maybe it could be interesting to mention some of the reference libraries that could be used in other languages to do bayesian inference (I guess stan - and its R and python interfaces - and pymc3 are possibly the most common libraries?).
| We need to choose prior probability distributions for each of the model parameters | ||
| that we are estimating. Let's use a normal distribution for the priors on | ||
| $\beta_{0}$ and $\beta_{1}$. The $\tau$ parameter must be a positive value | ||
| greater than 0 so let's use the gamma distribution as the prior on $\tau$. |
There was a problem hiding this comment.
maybe also a link to https://en.wikipedia.org/wiki/Gamma_distribution ?
| We cannot analytically solve the posterior probability distribution of our | ||
| linear model as the integration of the marginal likelihood is intractable. | ||
|
|
||
| But we can approximate it with markov chain monte carlo (MCMC) thanks to this |
There was a problem hiding this comment.
we could add also here a link to wikipedia of MCMC (or in the header below, maybe mentioning again the full meaning of the abbreviation).
| This is called the acceptance ratio. | ||
| 6. All proposals with greater probability than the current state are accepted. | ||
| So a proposeal is accepted if the acceptance ratio is greater than 1. If the | ||
| acceptance ratio is less than 1 then it is accepted with probability \alpha. In practice we can make a |
There was a problem hiding this comment.
| acceptance ratio is less than 1 then it is accepted with probability \alpha. In practice we can make a | |
| acceptance ratio is less than 1 then it is accepted with probability $\alpha$. In practice we can make a |
| nbImage("images/trace-sd.png") | ||
|
|
||
| nbText: hlMd""" | ||
| # Burnin |
There was a problem hiding this comment.
| # Burnin | |
| # Burn-in |
also below in various occurences. I looked for references about burn-in and found these, http://users.stat.umn.edu/~geyer/mcmc/burn.html https://www.johndcook.com/blog/2016/01/25/mcmc-burn-in/ not sure if there is something that could be said about it.
| nbText: hlMd""" | ||
| # Histograms | ||
| Now that we have our post burnin samples, let's see what our posterior probability | ||
| distributions look like for each model paramter. |
There was a problem hiding this comment.
| distributions look like for each model paramter. | |
| distributions look like for each model parameter. |
| geom_histogram(position="identity", alpha=some(0.5)) + | ||
| ggsave(fmt"{prefix}sd.png") | ||
|
|
||
| plotHistograms(@[burnin1, burnin2], "images/hist-") |
There was a problem hiding this comment.
what do the 1 and 2 represent in the plots? shouldn't we have a single distribution?
|
|
||
| nbText: hlMd""" | ||
| # Trace plots | ||
| We can get a sense for how well our mcmc performed and therefore gain some |
There was a problem hiding this comment.
| We can get a sense for how well our mcmc performed and therefore gain some | |
| We can get a sense for how well our MCMC performed and therefore gain some |
here and in other places occasionally mcmc is small caps (in code it has different meaning, in text it could be useful to refer to it always in the same way?)
| do this. | ||
|
|
||
| ## Combine Chains | ||
| Before we summarize the posterior. Let's combine the samples into one. |
There was a problem hiding this comment.
| Before we summarize the posterior. Let's combine the samples into one. | |
| Before we summarize the posterior, let's combine the samples into one. |
|
Ah maybe we could also change the title of the PR to Bayesian Inference or similar ;) (so it is easier to find in notifications) |
I finally got around to polishing up my tutorial on Bayesian Inference in Nim and would love to have it included in the SciNim project as I suggested in #19 (comment)