I need to fit a model using MLE with Frank copula by linking two discrete univate distribution function $u = F(x)$ and $v = F(y)$ together, and the joint distribution function is
$$ \Phi(x,y) = C(F(x),G(y)) $$
where
$$ C(u,v) = \frac{1}{\theta} \log\!\left( 1+\frac{(\exp(-\theta u)-1)(\exp(-\theta v)-1)}{\exp(-\theta)-1} \right) $$
and end up with the copula density as
$$\begin{align} c(F(x),G(y)) = &~C(F(x),G(y)) - C(F(x-1),G(y)) - C(F(x),G(y-1)) \\ & + C(F(x-1),G(y-1)) \end{align}$$
Therefore, I have to evaluate $\log c(u,v)$ and ultimately I have to evaluate $C$ and its first derivatives $C_u$, $C_v$ and $C_{\theta}$.
I couldn't get it work with my optimization algorithm (I am using fmin_bfgs and fmin_l_bfgs_b from scipy).
Most of the time, fmin_bfgs try to evaluate some very big point and make my code overflow (get -inf) and end with incorrect result.
I tried to use fmin_l_bfgs_b so that I can put bounds on the variable to prevent overflow, but I end up to only make it works (not overflow) in a too small interval, $[-2,2]$ or $[-1,1]$ in some cases. The resulting point will be on the bound and it must not be the correct solution.
I guess the problem is that $\log c$, $C_u$, $C_v$ and $C_{\theta}$ involve too much exp and log and I lost too much precision similar to the case like $\exp(x)-1$ and $\log(1+x)$.
Are there any better method to compute $C$, $C_u$, $C_v$ and $C_{\theta}$ accurately? (or compute $\log c$ and its derivatives directly?)
I am computing them using:
def frank_copula(u,v,t):
if t == 0.:
return u*v
else:
x1 = expm1(-u*t)
y1 = expm1(-v*t)
z1 = expm1(-t)
return log1p(x1*y1/z1)/(-t)
def frank_copula_u(u,v,t):
if t == 0.:
return v
else:
x = exp(t*u)
y = exp(t*v)
y1 = expm1(t*v)
z = exp(t)
return z*y1/(z*(x+y1) - x*y)
def frank_copula_v(u,v,t):
if t == 0.:
return u
else:
x = exp(t*u)
x1 = expm1(t*u)
y = exp(t*v)
z = exp(t)
return z*x1/(z*(x1+y) - x*y)
def frank_copula_t(u,v,t):
if t == 0.: # central difference here
h = 1e-8
return (frank_copula(u,v,t+h) - frank_copula(u,v,t-h))/(2.*h)
else:
x = exp(-t*u)
y = exp(-t*v)
z = exp(-t)
x1 = expm1(-t*u)
y1 = expm1(-t*v)
z1 = expm1(-t)
tmp1 = -z1*(x1*y*v + y1*x*u) + x1*y1*z
tmp2 = z1*z1 + x1*y1*z1
tmp3 = log1p(x1*y1/z1)/(-t)
return (tmp1/tmp2 + tmp3)/(-t)
EDIT
I fixed the dependency parameter $\theta$ first, and then find out the parameter of $F(x)$ and $G(y)$ together by maximizing $\log L(p) = \sum \log c_{\theta}(F(x|p),G(y|p))$. They involve the parameter $p$ of dimension 2 to 5.
Here is the output of how fmin_l_bfgs_b algorithm evaluate my (minus) log-likelihood function using my workaround, set the bounds to $[-0.2,1.2]$ and $[-1,1]$:
L: array([ 0., 0., 0., 0., 0.]) => 769.446206148
dL: array([ 0., 0., 0., 0., 0.]) => array([-338.93276474, -114.42079214, -216.55003298, -137.12475499,
-52.57668841])
L: array([ 1.2, 1.2, 1.2, 1. , 1. ]) => 8019.20848421
dL: array([ 1.2, 1.2, 1.2, 1. , 1. ]) => array([ 10020.77504911, 6095.3498353 , 8824.13266013, 6763.34035373,
2021.60815303])
L: array([ 0.31272257, 0.31272257, 0.31272257, 0.26060214, 0.26060214]) => 638.623452278
dL: array([ 0.31272257, 0.31272257, 0.31272257, 0.26060214, 0.26060214]) => array([-77.88688381, 54.50079579, 11.01938969, 98.73509615, 78.99317652])
L: array([ 0.41609651, 0.21822845, 0.28865935, 0.12751238, 0.13439021]) => 612.115679161
dL: array([ 0.41609651, 0.21822845, 0.28865935, 0.12751238, 0.13439021]) => array([-92.53682091, 12.87793807, -22.45693087, 26.74828638, 25.65511475])
L: array([ 7.05680743e-01, 1.28761254e-01, 3.21734086e-01,
-3.25757720e-02, 4.93345936e-04]) => 590.768526126
dL: array([ 7.05680743e-01, 1.28761254e-01, 3.21734086e-01,
-3.25757720e-02, 4.93345936e-04]) => array([-15.61044704, 16.3196277 , 16.76613757, 10.17553022, 1.73461386])
L: array([ 0.81933183, 0.05230577, 0.26232909, -0.07694669, -0.01632839]) => 587.63001378
dL: array([ 0.81933183, 0.05230577, 0.26232909, -0.07694669, -0.01632839]) => array([-10.36049322, 5.82433346, 10.12673749, -7.60549293, -3.70091331])
L: array([ 0.93612938, -0.03119036, 0.15964292, -0.05604783, -0.02092306]) => 585.655142497
dL: array([ 0.93612938, -0.03119036, 0.15964292, -0.05604783, -0.02092306]) => array([-1.60415082, 1.31457884, 5.95207277, -9.31825269, -5.19784526])
L: array([ 0.99490642, -0.09547364, 0.05431211, 0.01089551, 0.01278564]) => 584.829795091
dL: array([ 0.99490642, -0.09547364, 0.05431211, 0.01089551, 0.01278564]) => array([-0.44946283, -0.58748774, -1.01623299, -2.83092327, 1.18001099])
L: array([ 0.99380363, -0.10055615, 0.0527463 , 0.03080186, -0.00557434]) => 584.840769319
dL: array([ 0.99380363, -0.10055615, 0.0527463 , 0.03080186, -0.00557434]) => array([ 0.29741289, 0.67052685, -0.11044264, 1.88188906, -3.31504582])
L: array([ 0.99442739, -0.09768137, 0.05363196, 0.01954239, 0.00481047]) => 584.813957801
dL: array([ 0.99442739, -0.09768137, 0.05363196, 0.01954239, 0.00481047]) => array([-0.13879393, -0.04846326, -0.63399355, -0.79472897, -0.7858597 ])
L: array([ 0.99245896, -0.09809096, 0.05491007, 0.02184427, 0.00701989]) => 584.811245584
dL: array([ 0.99245896, -0.09809096, 0.05491007, 0.02184427, 0.00701989]) => array([ 0.16175504, 0.3075338 , -0.23277139, -0.09032161, -0.17666573])
L: array([ 0.99043781, -0.09901399, 0.056742 , 0.02308671, 0.00811308]) => 584.810330203
dL: array([ 0.99043781, -0.09901399, 0.056742 , 0.02308671, 0.00811308]) => array([ 0.08694667, 0.38532343, -0.11389679, 0.14966267, 0.06054983])
L: array([ 0.98934696, -0.10209532, 0.05898625, 0.02414812, 0.00863789]) => 584.809284488
dL: array([ 0.98934696, -0.10209532, 0.05898625, 0.02414812, 0.00863789]) => array([ 0.04661414, 0.30695628, -0.01908606, 0.20837613, 0.14123964])
L: array([ 0.9893386 , -0.10784629, 0.06164861, 0.02504006, 0.00849419]) => 584.808428091
dL: array([ 0.9893386 , -0.10784629, 0.06164861, 0.02504006, 0.00849419]) => array([ 0.00282033, 0.01810351, 0.00812252, 0.0266529 , 0.02348458])
L: array([ 0.98941539, -0.10801042, 0.06168044, 0.02499464, 0.00841819]) => 584.80842517
dL: array([ 0.98941539, -0.10801042, 0.06168044, 0.02499464, 0.00841819]) => array([-0.00429622, -0.00085152, -0.00194202, 0.00036292, 0.0010925 ])
and here is the output if I use my original bounds $[-5,5]$ (this bound works if I treat $F(x)$ and $G(y)$ as independent):
(Remarks: If I treat them as independent, I just need to evaluate the density function $f$,$g$ and they are much more stable and therefore I guess my problem is due to inaccuracy)
L: array([ 0., 0., 0., 0., 0.]) => 769.446206148
dL: array([ 0., 0., 0., 0., 0.]) => array([-338.93276474, -114.42079214, -216.55003298, -137.12475499,
-52.57668841])
L: array([ 5., 5., 5., 5., 5.]) => inf
dL: array([ 5., 5., 5., 5., 5.]) => array([ nan, nan, nan, nan, nan])
L: array([ 0., 0., 0., 0., 0.]) => 769.446206148
dL: array([ 0., 0., 0., 0., 0.]) => array([-338.93276474, -114.42079214, -216.55003298, -137.12475499,
-52.57668841])
and sorry I just found that fmin_bfgs works and will converge to the same point got by using my workaround method. But it will evaluate some big point and get -inf and nan:
dL: array([ 0., 0., 0., 0., 0.]) => array([-338.93276474, -114.42079214, -216.55003298, -137.12475499,
-52.57668841])
L: array([ 0., 0., 0., 0., 0.]) => 769.446206148
L: array([ 338.93276474, 114.42079214, 216.55003298, 137.12475499,
52.57668841]) => inf
dL: array([ 338.93276474, 114.42079214, 216.55003298, 137.12475499,
52.57668841]) => array([ nan, nan, nan, nan, nan])
L: array([ 3.38932765e-06, 1.14420792e-06, 2.16550033e-06,
1.37124755e-06, 5.25766884e-07]) => 769.444241862
dL: array([ 3.38932765e-06, 1.14420792e-06, 2.16550033e-06,
1.37124755e-06, 5.25766884e-07]) => array([-338.93116156, -114.42016384, -216.54896297, -137.12398571,
-52.57656892])
L: array([ 0.00238973, 0.00080675, 0.00152684, 0.00096683, 0.00037071]) => 768.063616272
dL: array([ 0.00238973, 0.00080675, 0.00152684, 0.00096683, 0.00037071]) => array([-337.80294205, -113.97709305, -215.79496551, -136.58082878,
-52.49109471])
L: array([ 169.46757724, 57.21079945, 108.27577991, 68.56286091,
26.28852956]) => inf
dL: array([ 169.46757724, 57.21079945, 108.27577991, 68.56286091,
26.28852956]) => array([ nan, nan, nan, nan, nan])
L: array([ 0.00238973, 0.00080675, 0.00152684, 0.00096683, 0.00037071]) => 768.063616272
dL: array([ 0.00238973, 0.00080675, 0.00152684, 0.00096683, 0.00037071]) => array([-337.80294205, -113.97709305, -215.79496551, -136.58082878,
-52.49109471])
L: array([ 338.93276474, 114.42079214, 216.55003298, 137.12475499,
52.57668841]) => inf
L: array([ 169.46638237, 57.21039607, 108.27501649, 68.56237749,
26.28834421]) => inf
L: array([ 84.73319119, 28.60519803, 54.13750824, 34.28118875, 13.1441721 ]) => inf
L: array([ 42.36659559, 14.30259902, 27.06875412, 17.14059437, 6.57208605]) => inf
L: array([ 21.1832978 , 7.15129951, 13.53437706, 8.57029719, 3.28604303]) => inf
L: array([ 10.5916489 , 3.57564975, 6.76718853, 4.28514859, 1.64302151]) => inf
L: array([ 5.29582445, 1.78782488, 3.38359427, 2.1425743 , 0.82151076]) => inf
L: array([ 2.64791222, 0.89391244, 1.69179713, 1.07128715, 0.41075538]) => inf
L: array([ 1.32395611, 0.44695622, 0.84589857, 0.53564357, 0.20537769]) => 2010.68946691
L: array([ 0.25288777, 0.08537274, 0.16157439, 0.10231284, 0.03922902]) => 651.514679659
dL: array([ 0.25288777, 0.08537274, 0.16157439, 0.10231284, 0.03922902]) => array([-214.23525183, -56.37608777, -123.86215026, -59.8381013 ,
-29.35530006])
L: array([ 2.45189744, -1.10216675, -0.00860584, -2.25562141, -0.0889141 ]) => 1140.89767258
dL: array([ 2.45189744, -1.10216675, -0.00860584, -2.25562141, -0.0889141 ]) => array([ 200.97116275, -107.57429828, 77.9175913 , -284.63382774,
222.78613298])
L: array([ 0.60625088, -0.10545527, 0.1342278 , -0.27658817, 0.01863746]) => 641.856514926
dL: array([ 0.60625088, -0.10545527, 0.1342278 , -0.27658817, 0.01863746]) => array([-179.74804571, -91.13312507, -130.64600443, -147.53451709,
-47.29257281])
L: array([ 0.94280274, -0.44136287, -0.34555214, -0.04312928, 0.3963997 ]) => 650.978041317
dL: array([ 0.94280274, -0.44136287, -0.34555214, -0.04312928, 0.3963997 ]) => array([-147.33106574, -91.10716372, -144.5277363 , -121.87325508,
59.49497911])
L: array([ 0.7272913 , -0.22626399, -0.03832441, -0.19262499, 0.15449914]) => 638.059222532
dL: array([ 0.7272913 , -0.22626399, -0.03832441, -0.19262499, 0.15449914]) => array([-170.21539216, -91.09157378, -134.59864658, -140.28045749,
-16.11667578])
...
L: array([ 0.98943719, -0.10800294, 0.06167076, 0.02498195, 0.00840848]) => 584.808425121
dL: array([ 0.98943719, -0.10800294, 0.06167076, 0.02498195, 0.00840848]) => array([ -6.05248630e-08, 2.30028962e-08, 3.00105580e-08,
-1.12310441e-07, 1.26875060e-08])
There may be a possibility that my gradient function is computed wrongly, therefore, I check the my gradient function againist finite difference:
def approx_gradient(f):
def wrapped_f(x0):
n = len(x0)
G = np.empty((n,),dtype=np.float64)
h = np.sqrt(np.finfo(np.float64).eps)
g1 = f(x0)
for i in xrange(n):
x1 = np.array(x0,dtype=np.float64)
x1[i] += h
g2 = f(x1)
G[i] = (g2 - g1) / h
return G
return wrapped_f
Print the error as:
# r1 is gradient computed by code, r2 is gradient computed by above finite-difference
d = np.log10(np.mean(np.abs(r1-r2)))
print("dL: %s => %s, error: %s\n" % (repr(x0),repr(r1),str(d)))
and the output:
dL: array([ 0., 0., 0., 0., 0.]) => array([-338.93276474, -114.42079214, -216.55003298, -137.12475499,
-52.57668841]), error: -1.12866497845
dL: array([ 338.93276474, 114.42079214, 216.55003298, 137.12475499,
52.57668841]) => array([ nan, nan, nan, nan, nan]), error: nan
dL: array([ 3.38932765e-06, 1.14420792e-06, 2.16550033e-06,
1.37124755e-06, 5.25766884e-07]) => array([-338.93116156, -114.42016384, -216.54896297, -137.12398571,
-52.57656892]), error: -1.3444523876
dL: array([ 0.00238973, 0.00080675, 0.00152684, 0.00096683, 0.00037071]) => array([-337.80294205, -113.97709305, -215.79496551, -136.58082878,
-52.49109471]), error: -1.76425541273
dL: array([ 169.46757724, 57.21079945, 108.27577991, 68.56286091,
26.28852956]) => array([ nan, nan, nan, nan, nan]), error: nan
dL: array([ 0.00238973, 0.00080675, 0.00152684, 0.00096683, 0.00037071]) => array([-337.80294205, -113.97709305, -215.79496551, -136.58082878,
-52.49109471]), error: -1.76425541273
dL: array([ 0.25288777, 0.08537274, 0.16157439, 0.10231284, 0.03922902]) => array([-214.23525183, -56.37608777, -123.86215026, -59.8381013 ,
-29.35530006]), error: -4.13744558124
dL: array([ 2.45189744, -1.10216675, -0.00860584, -2.25562141, -0.0889141 ]) => array([ 200.97116275, -107.57429828, 77.9175913 , -284.63382774,
222.78613298]), error: 4.55716763394
dL: array([ 0.60625088, -0.10545527, 0.1342278 , -0.27658817, 0.01863746]) => array([-179.74804571, -91.13312507, -130.64600443, -147.53451709,
-47.29257281]), error: -3.6510702955
dL: array([ 0.94280274, -0.44136287, -0.34555214, -0.04312928, 0.3963997 ]) => array([-147.33106574, -91.10716372, -144.5277363 , -121.87325508,
59.49497911]), error: -2.92020205735
dL: array([ 0.7272913 , -0.22626399, -0.03832441, -0.19262499, 0.15449914]) => array([-170.21539216, -91.09157378, -134.59864658, -140.28045749,
-16.11667578]), error: -3.52122115041
dL: array([ 0.87382717, -0.37020575, -0.21805711, -0.03796395, 0.05122782]) => array([-164.76068544, -94.42547784, -141.97757919, -126.19087232,
-39.8567946 ]), error: -3.19928906949
dL: array([ 0.87817125, -0.41141175, -0.15890429, -0.01557913, 0.05239629]) => array([-153.02101444, -89.44357123, -130.9150735 , -117.22771105,
-36.90142624]), error: -3.0131329785
dL: array([ 0.89554758, -0.57623572, 0.07770701, 0.07396015, 0.05707019]) => array([-98.94618778, -66.39110534, -79.97680149, -75.88852736, -21.65497929]), error: -3.5340970988
dL: array([ 1.31009037, 1.3391422 , -0.80615157, -0.3413817 , -0.09094827]) => array([ 88.75622027, 103.66551146, 28.52931607, 59.92216422,
34.82483184]), error: -3.86721832568
dL: array([ 1.03325948, 0.06005647, -0.2159125 , -0.0640172 , 0.00789817]) => array([-57.66593306, -27.77908949, -59.01722568, -47.64785975, -13.70147329]), error: -3.21832537521
dL: array([ 9.95932857e-01, -1.40930257e-01, 1.39440628e-01,
4.05132617e-02, 5.47288338e-04]) => array([ 27.98975787, 12.98469404, 25.68770109, 19.28943637, 5.48862252]), error: -4.19602531575
dL: array([ 0.98974304, -0.10584977, 0.05390259, 0.02224462, 0.00674051]) => array([-3.07471251, -1.62491769, -2.83340373, -2.41904387, -1.23502546]), error: -3.84209274725
dL: array([ 0.98982872, -0.10776008, 0.0609875 , 0.0246589 , 0.00848128]) => array([-0.07129277, -0.05980753, -0.10964847, -0.12157649, 0.00033618]), error: -3.57935797704
dL: array([ 0.98941444, -0.10804747, 0.06173506, 0.02499305, 0.00836133]) => array([-0.00140251, -0.00150752, 0.00292727, 0.00098813, -0.01230324]), error: -4.19996237842
dL: array([ 0.98944092, -0.10800269, 0.06166992, 0.0249806 , 0.00840967]) => array([ 0.00141414, 0.000485 , 0.00088848, 0.00042683, 0.00068224]), error: -4.05494671838
dL: array([ 0.98943748, -0.10800296, 0.06167061, 0.02498174, 0.00840839]) => array([ 2.32462795e-05, -1.47758964e-05, -5.04953544e-06,
-4.29179409e-05, -1.69174238e-05]), error: -4.30601047846
dL: array([ 0.9894372 , -0.10800295, 0.06167076, 0.02498195, 0.00840848]) => array([ 3.89230032e-07, -8.59830965e-07, -1.92379457e-07,
-1.03506122e-06, -1.88641715e-06]), error: -4.04379091415
dL: array([ 0.98943719, -0.10800294, 0.06167076, 0.02498195, 0.00840848]) => array([ -6.05248630e-08, 2.30028962e-08, 3.00105580e-08,
-1.12310441e-07, 1.26875060e-08]), error: -3.67971485116
OK. I think my updated question would be:
- I still think my gradient does not compute accurately in the box $[-5,5]^n$. I will approximate fisher information using my gradient function, I couldn't get a semi-positive definite information matrix.
- Are there any way to avoid BFGS to evaluate some "big" points?
- Are there any way to compute the copula density by computing $f$ and $g$ only for discrete margin?
