Problem with Ln?

[mikldk]

Consider (n) indpendent observations from exponential distribution with rate parameter (r). Let the sum of the observations be denoted by (S).

It is often interesting to work with the joint density function. First note that when (S) enters as being known, the joint density function is a function of the rate, (r), alone and this function is called the likelihood function. The maximum likelihood estimate of (r) is the value of (r) that maximizes the likelihood function (or equivalently, the log likelihood function).

In mathematics, letting (xi) for (i = 1, \ldots, n) be the observations and (S = \sum{i=1}^n x_i), then the log likelihood function becomes and the log likelihood function becomes

We want to work with this in yacas: take the log (as normally log likelihood functions are easier to work with as the product is turned into a sum) and find critical points by differentiating (where the sum is simpler to handle than the product).

## In> L := r^n * Exp(-r*S)
## Out> r^n*Exp(-r*S)

If we try this in yacas it cannot handle the fairly simple rules like (\log):

## In> Ln(L)
## Out> Ln(r^n*Exp(-r*S))

## In> Simplify(Ln(L))
## Out> Ln(r^n*Exp(-r*S))

## // Expected n*Ln(r) - r*S

And not in differentiation either:

## In> D(r) Ln(L)
## Out> (n*r^(n-1)*Exp(-r*S)-r^n*S*Exp(-r*S))/(r^n*Exp(-r*S))

## In> Simplify(D(r) Ln(L))
## Out> (n*r^(n-1)-r^n*S)/r^n

## // Expected n/r - S

It is only if we manually take the log we get the expected result:

## In> D(r) n*Ln(r)-r*S
## Out> n/r-S

[grzegorzmazur]

Hi,

As always, good catch :) Unfortunately, the whole thing is a bit more complicated, the complication stemming from the fact that there's no good (not to mention operational) universal notion of simple. Which form of an expression is better suited for particular calculations does depend on the context and their purpose.

Often the typical cases can be dealt with, but not always. This is the reason for which Simplify() does not touch your expression. But not all hope is lost :) There're two specialized function, LnExpand() and LnCombine(), with the first doing what you're after

In> LnExpand(Ln(L))
Out> Ln(r)*n-r*S

There are two issues though:

• the functions are not exported; this one is easy, I'll commit a fix later today
• application of LnExpand() often has to be followed by Simplify(), and the results may be only just passable

  In> LnExpand(D(r)Ln(L))
  Out> (n*r^(n-1)*Exp(-r*S)-r^n*S*Exp(-r*S))/(r^n*Exp(-r*S))
  In> Simplify(%)
  Out> (n*r^(n-1)-r^n*S)/r^n

Fortunately, at least in this case you can get better results doing

In> D(r)LnExpand(Ln(L))
Out> n/r-S

For now I'll export the LnExport() and LnCombine() functions; improving Simplify() is on my todo list, but it's not going to be an easy task and I can't promise any timeline.

Anyway, thanks a lot for all the work you do, getting real-world examples of problems with yacas helps me immensely with improving it.

Cheers, Grzesiek

[mikldk]

Thank you for the very nice explanation (and I if course agree with your "simple" discussion). No worries about the timeline.

Cheers, Mikkel.