3

I am trying to implement an empirical bayesian ML-II（Maximum likelihood estimation Type II）method for estimating prior distribution parameters from historical data

Where：

  1. π(θ) is an expression for the prior distribution
  2. p(x|θ) is is an expression for the data distribution
  3. m(x) is an expression for the marginal distribution

According to the steps, I need to first integrate to find the expression of the marginal distribution, and then find the extreme value of this expression to estimate the parameters of the prior distribution. Extreme values can be achieved using methods such as scipy.optimize. So the question is how do we integrate this?

enter image description here

abraxas
  • 45
  • 9
  • @merv Our Dear Merv,Thank you. I know of packages for bayesian models like pymc3.I asked them, and they said that empirical bayes is not a universal method, and they didn't implement ML-II.In particular, this ML-II is not a MAP.It looks like the pymc3 method for MAP doesn't do ML-II – abraxas Mar 04 '19 at 14:41
  • Perhaps `symfit` will be able to help you, have a look at [this](https://symfit.readthedocs.io/en/stable/fitting_types.html#likelihood) example in the docs. Then you will be able to use `sympy` style analytical expressions, and fit them to your data using `scipy` but without having to interact with `scipy`. Disclaimer: I'm the author of `symfit`. – tBuLi Mar 04 '19 at 15:43
  • @tBuLi I have read your link, it seems that I still need to write the formula of model by myself, and this is exactly the problem I encountered, I need to integrate to find the expression of marginal distribution first.If I were to use my brain integral to figure out this expression, which is where I'd get into trouble, I'd like to be able to write some generic code to integrate this expression – abraxas Mar 05 '19 at 06:39
  • Am I correct in thinking that you have an analytical expression for your model, but would like to have some automated way to integrate it into a marginal distribution, which you then need to fit to a dataset? Let me know if that is the case, then I will write you an example of how to do that with symfit. – tBuLi Mar 05 '19 at 17:27
  • @tBuLi I am impressed by your understanding! Indeed, the method to estimate the parameters of the prior distribution is to integrate the prior distribution multiplied by the likelihood function, and then find out what kind of prior distribution parameters can maximize the value of the marginal distribution after obtaining the expression of the marginal distribution.Then, the parameters will be used as the parameters of the prior distribution to carry out the maximum posterior estimation – abraxas Mar 06 '19 at 03:02

1 Answers1

2

Here's a go at this using symfit. As an example I choose to sample from a bivariate normal distribution with no covariance.

import numpy as np
import matplotlib.pyplot as plt
from symfit import Model, Fit, Parameter, Variable, integrate, oo
from symfit.distributions import Gaussian
from symfit.core.objectives import LogLikelihood

# Make variables and parameters
x = Variable('x')
y = Variable('y')
m = Variable('m')
x0 = Parameter('x0', value=0.6, min=0.5, max=0.7)
sig_x = Parameter('sig_x', value=0.1)
y0 = Parameter('y0', value=0.7, min=0.6, max=0.9)
sig_y = Parameter('sig_y', value=0.05)

pdf = Gaussian(x=x, mu=x0, sig=sig_x) * Gaussian(x=y, mu=y0, sig=sig_y)
marginal = integrate(pdf, (y, -oo, oo), conds='none')
print(pdf)
print(marginal)

model = Model({m: marginal})

# Draw 10000 samples from a bivariate distribution
mean = [0.59, 0.8]
cov = [[0.11**2, 0], [0, 0.23**2]]
xdata, ydata = np.random.multivariate_normal(mean, cov, 10000).T

# We provide only xdata to the model
fit = Fit(model, xdata, objective=LogLikelihood)
fit_result = fit.execute()
print(fit_result)

xaxis = np.linspace(0, 1.0)
plt.hist(xdata, bins=100, density=True)
plt.plot(xaxis, model(x=xaxis, **fit_result.params).m)
plt.show()

This prints the following for the pdf and the marginal distribution:

>>> exp(-(-x0 + x)**2/(2*sig_x**2))*exp(-(-y0 + y)**2/(2*sig_y**2))/(2*pi*Abs(sig_x)*Abs(sig_y))
>>> sqrt(2)*sig_y*exp(-(-x0 + x)**2/(2*sig_x**2))/(2*sqrt(pi)*Abs(sig_x)*Abs(sig_y))

And for the fit results:

Parameter Value        Standard Deviation
sig_x     1.089585e-01 7.704533e-04
sig_y     5.000000e-02 nan
x0        5.905688e-01 -0.000000e+00
Fitting status message: b'CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH'
Number of iterations:   9
Regression Coefficient: nan

enter image description here

You can see that x0 and sig_x have been obtained properly, but no information can be obtained about the parameter to do with y. I think that makes sense in this example since there is no correlation, but I'll leave you to struggle with those kinds of details ;).

tBuLi
  • 2,295
  • 2
  • 16
  • 16
  • Our Dear Li.You left a deep impression on me!But there are some differences here: if y0 is a data distribution, then its mean is subject to a gaussian distribution. It seems that the mean of y0 is a fixed parameter, not a random variable subject to a gaussian distribution. How can I modify your code to achieve this? – abraxas Mar 07 '19 at 05:16
  • Stop it, you are making me blush. Yes, you are right, in my example I did used two uncorrelated Gaussians, so `y0` has no influence on the fit since it does not appear in the marginal distribution. But if they are dependent on each other, `y0` it will not be integrated away. So in order to do your model, you have to adapt my `pdf` to reflect that for your problem. So: `pdf = pi * p`, and then integrate out `theta`. That should give you `m`. – tBuLi Mar 07 '19 at 14:23