Message Boards Message Boards

How do I solve a profile interval (statistics)

Posted 10 years ago

I want to calculate a 95% profile interval for f0. Do worry too much if you do not understand what it is. I will illustrate with a simple example.

mylik[f0_,p_]:=Module[
    {ans},
    ans = 
    - Log[2] 
    - Log[120] 
    - Log[6227020800] 
    - Log[1124000727777607680000] 
    - Log[15511210043330985984000000] 
    + f0 Log[(1 - p)^6] 
    + 25 Log[6 (1 - p)^5 p] 
    + 22 Log[15 (1 - p)^4 p^2] 
    + 13 Log[20 (1 - p)^3 p^3] 
    + 5 Log[15 (1 - p)^2 p^4] 
    + Log[6 (1 - p) p^5] 
    + 2 Log[p^6] 
    - LogGamma[1 + f0] 
    + LogGamma[69 + f0]
];
{maxlog, mles} = NMaximize[{mylik[f0, p], f0 > 0, 0 < p < 1}, {f0, p}];
proflikfun[f0_] := Module[
    {ans, p},

    ans = NMaximize[{mylik[f0, p], 0 < p < 1}, p];
    ans = ans[[1]];
    ans = ans + Quantile[ChiSquareDistribution[1], 0.95]/2 - maxlog;
    ans
];

Now to plot the function. Because I know the answer, so I will plot the range

Plot[proflikfun[f0], {f0, 0, 20}];

enter image description here

Now, how do I solve this function at 0?

I tried things like

FindRoot[proflikfun[f0] == 0, {f0, 0}]
NSolve[proflikfun[f0] == 0, {f0, 0}]
POSTED BY: Casper YC
4 Replies

The error message you see says that something is being passed a symbol instead of a number and that it needs a number instead of a symbol. Sometimes it is important to make sure that functions only evaluate when you give them a number and won't evaluate when given them a symbol. You can have functions only evaluate when given a number by using a pattern,?NumericQ. See this piece of documentation here http://support.wolfram.com/kb/12502.

mylik[f0_?NumericQ, p_?NumericQ] := 
  Module[{ans}, 
   ans = -Log[2] - Log[120] - Log[6227020800] - 
     Log[1124000727777607680000] - Log[15511210043330985984000000] + 
     f0 Log[(1 - p)^6] + 25 Log[6 (1 - p)^5 p] + 
     22 Log[15 (1 - p)^4 p^2] + 13 Log[20 (1 - p)^3 p^3] + 
     5 Log[15 (1 - p)^2 p^4] + Log[6 (1 - p) p^5] + 2 Log[p^6] - 
     LogGamma[1 + f0] + LogGamma[69 + f0]];
{maxlog, mles} = NMaximize[{mylik[f0, p], f0 > 0, 0 < p < 1}, {f0, p}];

proflikfun[f0_?NumericQ] := 
  Module[{ans, p}, ans = NMaximize[{mylik[f0, p], 0 < p < 1}, p];
   ans = ans[[1]];
   ans = ans + Quantile[ChiSquareDistribution[1], 0.95]/2 - maxlog;
   ans];

FindRoot uses normally uses things like Newton's method. It differentiates its input which can't be done in this case. It can however use secant like methods which don't use derivatives if you give it two starting points.

FindRoot[proflikfun[f0], {f0, 1, 2}]

NSolve uses a variety of methods but doesn't take a starting point. Instead, you list what variables to solve for. The correct syntax is:

NSolve[proflikfun[f0] == 0, f0]
POSTED BY: Sean Clarke
Posted 10 years ago

That has been really helpful. I still wonder, is it possible to get both roots (if the function has two)? In this case, I would expect two, but both NSolve and FindRoot only give me one.

I looke at the help page on FindRoot, which says, " .......searches for a numerical root of f........", so FindRoot is out of the question?

It is possible to use NSolve?

POSTED BY: Casper YC
POSTED BY: Sean Clarke
Posted 10 years ago

Thanks Sean again for the detailed explanation.

Unfortunately, this is only a 'simpler' example to illustrate what originally i was struggling. My real problems deals with bigger models with lots of other parameters, not just two.

Note that in the most recent reply, ?NumericQ was note used. That was part of the reason why I did not use it in the first instance. Apart from the likelihood estimates, I was also trying to use D[] to work out hessian (2nd derivative), and ect.

mylik[f0_?NumericQ, p_?NumericQ] := -Log[2] - Log[120] - 
  Log[6227020800] - Log[1124000727777607680000] - 
  Log[15511210043330985984000000] + f0 Log[(1 - p)^6] + 
  25 Log[6 (1 - p)^5 p] + 22 Log[15 (1 - p)^4 p^2] + 
  13 Log[20 (1 - p)^3 p^3] + 5 Log[15 (1 - p)^2 p^4] + 
  Log[6 (1 - p) p^5] + 2 Log[p^6] - LogGamma[1 + f0] + 
  LogGamma[69 + f0]
FullSimplify[D[mylik[f0, p], p] == 0]

(will not work with ?NumericQ )

In order to build up a larger function (programm), i did not use ?NumericQ.

But again, your input is really helpful to me as I have only been using MMA substantially for just about 3 weeks.

For now,

FindRoot[proflikfun[f0], {f0, 0, 10}]
FindRoot[proflikfun[f0], {f0, 10, 20}]

works fine. I will try to get around with more lines, and hopefully generalize it for my other data sets and models.

POSTED BY: Casper YC
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract