On Jan 2, 2015, at 12:04 PM, Ben Bolker wrote:
> Ben Bolker <bbolker <at> gmail.com> writes:
>
>>
>> Paul Hudson <paulhudson028 <at> gmail.com> writes:
>>
>
> [snip]
>
>> library("tweedie")
>> set.seed(1001)
>> r <- rtweedie(1000,1.5,mu=2,phi=2)
>> library("bbmle")
>> dtweedie2 <- function(x,power,mu,phi,log=FALSE,debug=FALSE) {
>> if (debug) cat(power,mu,phi,"\n")
>> res <- dtweedie(y=x,xi=power,mu=mu,phi=phi)
>> if (log) log(res) else res
>> }
>> m <- mle2(r~dtweedie2(power=exp(logpower),
>> mu=exp(logmu),
>> phi=exp(logphi)),
>> ## don't start with logpower=0 (power=1)
>> start=list(logpower=0.1,logmu=0,logphi=0),
>> data=data.frame(r),
>> method="Nelder-Mead")
>>
>> dtweedie2(r,xi=exp(0.1),mu=1,phi=1)
That threw an error.
>>
>> In principle MASS::fitdistr could be made to work too.
>
> PS in hindsight, you're better off with the built-in tweedie.power()
> recommended by another poster. Estimating the power parameter for
> Tweedie distributions is known to be difficult, and the naive approach I show
> above may only work in best-case scenarios.
I did try fitdistr in the ftdistrplus package multiple errors due to
estiamating power parameters less than 1.00 before I noticed the
tweedie.profile function. (Your implementation via a transformation effectively
sidestepped that possibility.) And then I kicked myself for not reading the
documentation more thoroughly.
Comparing the two methods:
> exp( m@coef )
logpower logmu logphi
1.516015 1.935660 2.089286
> tweedie.profile(r~1)[c('xi.max', 'phi.max')]
1.2 1.3 1.4 1.5 1.6 1.7 1.8
.......Done.
$xi.max
[1] 1.514286
$phi.max
[1] 2.085714
The bbmle::mle2 approach is about 25 times faster than 'tweedie.profile'.
--
David Winsemius
Alameda, CA, USA
______________________________________________
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
