I am exploring how to model a data set using normal distributions with both mean and variance defined as linear functions of independent variables.
Something like N ~ (f(x), g(x)).
I generate a random sample like this:
def draw(x): return norm(5 * x + 2, 3 *x + 4).rvs(1)
So I want to retrieve 5, 2 and 4 as the parameters for my distribution.
I generate my sample:
smp = np.zeros((100,2))
for i in range(0, len(smp)): smp[i] = i smp[i] = draw(i)
The likelihood function is:
def lh(p): p_loc_b0 = p p_loc_b1 = p p_scl_b0 = p p_scl_b1 = p l = 1 for i in range(0, len(smp)): x = smp[i] y = smp[i] l = l * norm(p_loc_b0 + p_loc_b1 * x, p_scl_b0 + p_scl_b1 * x).pdf(y) return -l
So the parameters for the linear functions used in the model are given in the p 4-variable vector.
Using scipy.optimize, I can solve for the MLE parameters using an extremely low xtol, and already giving the solution as the starting point:
fmin(lh, x0=[2,5,3,4], xtol=1e-35)
Which does not work to well:
Warning: Maximum number of function evaluations has been exceeded. array([ 3.27491346, 4.69237042, 5.70317719, 3.30395462])
Raising the xtol to higher values does no good.
So i try using a starting solution far from the real solution:
>>> fmin(lh, x0=[1,1,1,1], xtol=1e-8) Optimization terminated successfully. Current function value: -0.000000 Iterations: 24 Function evaluations: 143 array([ 1., 1., 1., 1.])
Which makes me think:
PDF are largely clustered around the mean, and have very low gradients only a few standard deviations away from the mean, which must be not too good for numerical methods.
So how does one go about doing these kind of numerical estimation in functions where gradient is very near to zero away from the solution?