0
|
4604 Views
|
4 Replies
|
0 Total Likes
View groups...
Share

# How do I solve a profile interval (statistics)

Posted 9 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 - Log - Log - Log - Log + 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[]; ans = ans + Quantile[ChiSquareDistribution, 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}]; 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}] 
4 Replies
Sort By:
Posted 9 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 - Log - Log - Log - Log + 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 9 years ago
 These functions are numeric solvers. They don't have any way of knowing about the properties of the function proflikfun. So they don't know how many roots it has or even if the function is well behaved. NSolve will try it's best to find all of the roots, but there's guarantee it can. You can always use FindRoot with many different starting values and see if any of them converge. There's no guarantees. Another method might be to rephrase it as a minimization problem and use a global minimizer like NMiniIt may be possible to derive a symbolic version of proflikfun. If it's a polynomial or something simple, then there might be a symbolic method to find all of its roots and then Solve might work. Mylik doesn't need to be defined with ?NumericQ or Module. The problem seems to simplify down pretty substantially: mylik[f0_, p_] := -Log - Log - Log - Log - Log + 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] For a fixed value of f0, when is mylik maximized? FullSimplify[D[mylik[f0, p], p] == 0] For a given value of f0, which value of p maximizes mylik? Solve[FullSimplify[D[mylik[f0, p], p] == 0], p] {{p -> 145/(6 (68 + f0))}} So we can get rid of a lot of the numerical algorithms in this code. That might be helpful.
Posted 9 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 9 years ago
 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 - Log - Log - Log - Log + 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[]; ans = ans + Quantile[ChiSquareDistribution, 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]