# FindMinimum returns different values for 2-variable symmetric log-likelihood function

I have data (purchase data, currently, one user, one item) that I generate according to a probability function that is a product of two logistic functions, i.e. $p = \frac{e^{b1}}{1+e^{b1}} * \frac{e^{b2}}{1+e^{b2}}$

Using the data generated I know how many purchases there were and then I write my likelihood function $p^k *(1-p)^{n-k}$ where k is the number of purchases of a certain item and n is the number of trials.

I then use Mathematica's FindMinimum on the log-likelihood function, which is simply -$Log[p^k *(1-p)^{n-k}]$ to find the minimum. This is analogous to calculating the MLE, which tells us that the minimum is achieved at the correct value p, which in our case would be the correct values b1,b2.

Notice the function is symmetric in b1 and b2, so we would expect the results of FindMinimum to be symmetric as well. Not only are they not always symmetric, they are not the values I used to generate the data, which is what I would expect to get.

An example run:

In:= FindMinimum[-Log[(Exp[b1]*
Exp[b2]/((1 + Exp[b1]) (1 + Exp[b2])))^3829 * (1 -
Exp[b1]*Exp[b2]/((1 + Exp[b1]) (1 + Exp[b2])))^6171], {b1,
0.7}, {b2, 0.89}]
Out:= {6654.66, {b1 -> 0.373932, b2 -> 0.603008}}


And when I use the values I expect: b1=b2=1, I get a symmetric answer, but it's nowhere near what I expect:

{6654.66, {b1 -> 0.484412, b2 -> 0.484412}}


I ran the above line of code with various values of k, generated by repeating the experiment and the same problems persist.

Is there something really obvious that I am missing? Can anyone help me sort this out? Thank you very much in advance!

• What were the parameter values that you used to creat the data?
– gwr
Mar 3 '17 at 10:15
• @gwr I usedb1=1, b2=1 Mar 3 '17 at 10:56
• Essentially you are running 10 000 independent Bernoulli-Trials given $p$ and I have generated lots of samples of that size using your formula for $p$. I don't come anywhere near the values for purchases you have gotten, In my experiments all works out as should and (1,1) is on a "ridge" of most likely value combinations...
– gwr
Mar 3 '17 at 11:02
• "Is there something really obvious that I am missing?" Yes. $p$ can be estimated but there is no data collected that you've described that allows you to separately estimate $b_1$ and $b_2$. Same thing if $p=a_1 a_2$ or $p=a b^c+\cos(d)$. Now if you had some covariates in the logistic equations that varied appropriately from trial to trial, then you could probably estimate the coefficients that define $p$.
– JimB
Mar 3 '17 at 14:18
• @gwr I think one's prior on the given functional form would have to be awfully good, too. But you raise an interesting point. A few simulations would be interesting in a theoretical way. However, if this problem was about weapons or some pharmaceutical product, I'd hate to be on the receiving end of the results with such tenuous assumptions.
– JimB
Mar 3 '17 at 17:04

To me it looks like there may be something wrong with your data generating process.

Essentially you are doing 10,000 trials where the result for each trial $X_i$ is distributed i.i.d. Bernoulli:

$\begin{equation} X_i \sim Bernoulli\big( p = \text{sig}(b1) \cdot \text{sig}(b2) \big) \end{equation}$

where $\text{sig}$ is the sigmoid function. This means that the number of purchases $Y|n$ given a sample size of $n$ purchases will follow a Binomial Distribution:

$Y|n=10,000 \sim Binomial( n, p )$

Thus

dist = With[ {b1 = 1, b2 = 1, n = 10000} ,
BinomialDistribution[ n , 1 / ( (1 + Exp[-b1]) (1 + Exp[-b2]) ) ]
]; (* this really saves the work of going through 10000 trials... *)


Now we can do this experiment itself lots of times and at the same time find the most likely parameter value-combinations:

(* we now do 10 x 10.000 trials thus 10 x a Binomial-experiment *)
(* for each experiment we look at the MLE for (b1,b2) but we are finding *)
(* combinations by letting b1 be fixed over some range and minimizing for b2 *)

sols = With[ { n = 10000, data = RandomVariate[ dist, 10 ] },
Map[

(* for each experiment we table combinations *)

ParallelTable[

(* we are only interested in the parameter values thus NArgMin *)

NArgMin[
{
-Log[
(Exp[b1]*Exp[b2]/((1 + Exp[b1]) (1 + Exp[b2])))^# *
(1 - Exp[b1]*Exp[b2]/((1 + Exp[b1]) (1 + Exp[b2])))^(n - #)
],
b1 == β1 && 0 <= b2 (* fixing b1 *)
},
{b1, b2},
Method -> "NelderMead" (* making this explicit *)
],
{ β1, Range[ 0, 1.2, 0.01 ] } (* fixed range for b1 *)
]&,
data
]
];

ListLinePlot[ sols,
AxesLabel -> { "b1", "b2" },
Epilog -> { Red, PointSize[ Large ], Point[ {1,1} ] }
] As you can see the "true" parameter values are on the line that contains all value-combinations.

We may also check this more formally by asking which parameter combinations will give the same probabilty for a purchase:

Reduce[
1/((1 + Exp[-b1])*(1 + Exp[-b2])) == p, {b1, b2},
Reals
] // FullSimplify


0 < p < 1 && b1 > Log[-(p/(-1 + p))] && b2 == Log[-(((1 + E^b1) p)/(E^b1 (-1 + p) + p))]

If we make a table of these parameter combinations and plot them, we see, that the shape of that iso-probability-curve matches the empirical ones obtained above rather closely:

combinations = Rest @ Quiet @ With[
{ p = E^2 / (1+E)^2 }, (* that is your true combination *)
Table[
{ b1, Log[ -p (1 + Exp[b1])/( p + (p - 1) Exp[b1] )] },
{ b1, Log[ -p / (p-1)], 1.2, 0.001 }
]
];

ListLinePlot[ combinations, PlotRange -> { {0, 1.2}, {0,7} } ] • Sorry, I just saw this. First of all, thank you very much, you are taking a lot of time and effort to help me out and I really appreciate it. Secondly, I don't follow... Mainly, the syntax, are you writing in Mathematica? And I'm not clear on what your conclusion is... I don't write fluent Mathematica so it will take me a while to decode your code :) but I will do it now. Mar 3 '17 at 12:28
• @Sareet Whatever is marked as code (gray) is Mathematica. Your problem essentially starts with your data generating process as the probability of accepting a purchase will be $1-p$ and not $p$ if you use uniform distribution for each trial and accept a purchase if the random number is higher than the given probability.
– gwr
Mar 3 '17 at 12:32
• Thank you, I am still trying to figure out exactly what you did, I'm on my way :) I generate my data with r in U(0,1) and then if r<=p I accept the purchase, if I understand you correctly then this corresponds to p^#*p^(n-#). Moreover, it seems all of this is done behind the scenes by Mathematica Mar 3 '17 at 12:41
• @grw, I think I understand - you're saying that there are infinitely many pairs of (b1, b2) that satisfy $sig(b_1)*sig(b_2) = p$. Mar 3 '17 at 13:16
• @Sareet Yes, exactly and that is why - depending on where you start searching - with numerical optimization you can end up somewhere on that line.
– gwr
Mar 3 '17 at 13:20

"Is there something really obvious that I am missing?" Yes.

p can be estimated but there is no data collected that you've described that allows you to separately estimate $b_1$ and $b_2$. Same thing if $p=a_1 a_2$ or $p=a b c + \cos(d)$. Now if you had some covariates in the logistic equations that varied appropriately from trial to trial, then you could probably estimate the coefficients that define $p$.

For example, suppose you had covariates $u_i$ and $v_i$ for the $i$-th observation and $p_i$ was defined as follows: Then $b_1$ and $b_2$ could be estimated.

If you have no covariates such that $p$ potentially changes for each observation, then you can just estimate $p$. The point is that the data that you describe (0 or 1 for each observation) does not allow you to estimate coefficients from functions like what you have or even to suggest what the underlying function is.

Update: Why you can't get there from here

The equivalent of finding the parameter values that maximize the likelihood (or minimize minus the likelihood) in this case is to find where the derivative of the log likelihood is zero. You have the log of the likelihood being

logL = k Log[p] + (n - k) Log[1 - p]


To find the value of p to maximize this function we can use the following:

Solve[D[logL, p] == 0, p]
(* {{p -> k/n}} *)


We see that given the only data you have ($k$ = the number of purchases of a certain item and $n$ is the number of trials) there is an explicit solution for p.

That's all the information you have. One can't be fanciful and then assume $p=a*b$ or $p$ is some more complex function and be able to estimate the parameters in that function. It's not a matter as to if k and n are large enough. There just isn't the kind of information collected to allow a more in-depth assessment of what $p$ is made of.

Maybe here's an appropriate analogy: Suppose you have two very large bags of coins. One bag is filled with coins each having the same probability $p$ of coming up heads. The other bag is filled with coins that all vary in the probability of heads but the average probability is $p$.

If we sample $n$ coins from one bag and toss them one at a time and record the total number of heads ($k$), that information is insufficient to be able to tell from which bag we are sampling.

The bags are certainly different. But with the information as described above we simply can't tell the difference. There are ways to sample from a bag such that we can infer from which bag we are sampling. (For example, as we choose a coin we could toss that coin multiple times and record the number of heads for each coin.)

To repeat ad nauseam: you can't get there from here with just $k$ and $n$.

• There are indeed good reasons for "Occam's Razor" and the call for parsimony I guess.
– gwr
Mar 3 '17 at 17:23
• Thank you very much for your answer. There are no covariates, this is a simple experiment with all parameters independent of each other and given the purchase data I am interested in finding the probablities of purchase for each item. I will start out with a seed value of $\frac{# purchases}{number of experiments}$ from which I can choose an initial (b1, b2) pair of values and then I need to think of a criterion for improving the expression for the probability (since clearly the previous fraction is not an accurate expression) Mar 5 '17 at 8:03