PEtab-dev / petab_test_suite

PEtab test suite
BSD 3-Clause "New" or "Revised" License
6 stars 2 forks source link

challenge for 0007 #11

Closed marcusrosenblatt closed 4 years ago

marcusrosenblatt commented 4 years ago

Hi,

for test model 0007, I get the correct chi2 value, but the likelihood differs.

As you can see in this R code snippet:

chi2 <- obj(pouter)$value chi2 [1] 3.226562 -0.5(chi2+sum(log(2pi*c(0.5,0.6,0.7)^2))) [1] -2.809449

dMod obtains -2.8 instead of -1.38 as in the given solution.

All other test cases obtain the correct results for both chi2 and llh. Presumably it has something to do with the different log-transformations. Could it be that you internally rescale the noise parameter for log10 or log. But I think this is wrong.

Best, Marcus

marcusrosenblatt commented 4 years ago

Hi,

ok, this case is now splitted into two different ones. But the inconsistency in the calculation of the likelihood is still there.

Let us proceed with the discussion for example 0007. Here we now have two observations obs_a lin 0.5 and obs_b log10 0.6.

chi2 is equal to 0.2683. The "true" solution of the likelihood is -0.54491.

Even with rescaling of the noise parameter I, however, get a different solution. Assuming noise parameters have to be rescaled, 0.6 is the noise of the nominal value (0.5714281). Rescaling gives an error of 0.6/(0.57 * log(10))=0.457. So we should have for the likelihood:

llh = -0.5(chi2 + log(2 pi 0.5^2) + log (2 pi * 0.457^2)) = -0.496.

Where is my mistake?

yannikschaelte commented 4 years ago

Hi Marcus, sorry I missed this issue.

what we assume for measurement b in 0007 (0016) is that the noisy data are log(10)-normally distributed https://en.wikipedia.org/wiki/Log-normal_distribution. Thus, we assume that if we take the (10er) logarithm of simulation and measurement, then that is a normal distribution with sigma equal the (not-scaled!) noise parameter.

What we effectively compute for this situation here is for a

0.5*((0.2-0.429)**2 / 0.5**2 + np.log(2*np.pi*0.5**2))

and for b

0.5*((np.log10(0.8)-np.log10(0.571))**2/0.6**2 + np.log(2*np.pi*0.6**2*0.8**2))

together this gives the 0.54491. Note in particular the additional (noise=0.8)**2 in the log term for b. This comes from the transformation rule, to ensure that the density of y is a proper random variable (see the densitiy in the wikipedia article). Might this not be in your calculation yet?

yannikschaelte commented 4 years ago

More specifically, you find here https://github.com/PEtab-dev/PEtab/blob/develop/petab/calculate.py#L356 the in total six different noise models.

yannikschaelte commented 4 years ago

On second thoughts (thanks for raising this point actually ...) I think we have an error in calculating the density in AMICI ourselves. The log-term should not be np.log(2*np.pi*0.6**2*0.8**2)) but np.log(2*np.pi*0.6**2*0.8**2*np.log(10)**2)). This will be updated soon.