[R] Using betareg package to fit beta mixture with given initial parameters
Achim Zeileis
Achim.Zeileis at uibk.ac.at
Wed Mar 22 11:34:01 CET 2017
On Wed, 22 Mar 2017, Michael Dayan wrote:
> The method of setting the initial values given lambda, alpha1, etc. should
> not depend on the exact values of lambda, alpha1, etc. in my situation,
> i.e. it does not depend on my data.
Presently, flexmix() that betamix() is built on cannot take the parameters
directly for initialization. However, it is possible to pass a matrix with
initial 'cluster' probabilities. This can be easily generated using
dbeta().
For a concrete example consider the data generated in this discussion on
SO:
http://stats.stackexchange.com/questions/114959/mixture-of-beta-distributions-full-example
Using that data with random starting values requires 42 iterations until
convergence:
set.seed(0)
m1 <- betamix(y ~ 1 | 1, data = d, k = 2)
m1
## Call:
## betamix(formula = y ~ 1 | 1, data = d, k = 2)
##
## Cluster sizes:
## 1 2
## 50 100
##
## convergence after 42 iterations
Instead we could initialize with the posterior probabilities obtained at
the observed data (d$y), the true alpha/beta parameters (10; 30 vs. 30;
10) and the true cluster proportions (2/3 vs. 1/3):
p <- cbind(2/3 * dbeta(d$y, 10, 30), 1/3 * dbeta(d$y, 30, 10))
p <- p/rowSums(p)
This converges after only 2 iterations:
set.seed(0)
m2 <- betamix(y ~ 1 | 1, data = d, k = 2, cluster = p)
m2
## Call:
## betamix(formula = y ~ 1 | 1, data = d, k = 2, cluster = p)
##
## Cluster sizes:
## 1 2
## 100 50
##
## convergence after 2 iterations
Up to label switching and small numerical differences, the parameter
estimates agree. (Of course, these are on the mu/phi scale and not
alpha/beta as explained in the SO post linked above.)
coef(m1)
## (Intercept) (phi)_(Intercept)
## Comp.1 1.196286 3.867808
## Comp.2 -1.096487 3.898976
coef(m2)
## (Intercept) (phi)_(Intercept)
## Comp.1 -1.096487 3.898976
## Comp.2 1.196286 3.867808
> On Mar 22, 2017 04:30, "David Winsemius" <dwinsemius at comcast.net> wrote:
>
>
>> On Mar 21, 2017, at 5:04 AM, Michael Dayan <mdayan.research at gmail.com>
> wrote:
>>
>> Hi,
>>
>> I would like to fit a mixture of two beta distributions with parameters
>> (alpha1, beta1) for the first component, (alpha2, beta2) for the second
>> component, and lambda for the mixing parameter. I also would like to set a
>> maximum of 200 iterations and a tolerance of 1e-08.
>>
>> My question is: how to use the betareg package to run the fit with initial
>> values for the parameters alpha1, beta1, alpha2, beta2 and lambda? I saw
> in
>> the documentation that I would need to use the 'start' option of the
>> betareg function, with start described as "an optional vector with
> starting
>> values for all parameters (including phi)". However I could not find how
> to
>> define this list given my alpha1, beta1, alpha2, beta2 and lambda.
>>
>> The current code I have is:
>> mydata$y <- <my_data>
>> bmix <- betamix(y ~ 1 | 1, data = mydata, k = 2, fstol = 1e-08, maxit =
>> 200, fsmaxit = 200)
>>
>>
>> And I suspect I would need to do something along the lines:
>>
>> initial.vals <- c(?, ?, ?, ?, ?)
>> bmix <- betamix(y ~ 1 | 1, data = mydata, k = 2, fstol = 1e-08, maxit =
>> 200, fsmaxit = 200, control=betareg.control(start=initial.vals)))
>>
>> But I do not know what to use for initial.vals.
>
> If there were sensitivity to data, then wouldn't that depend on your
> (unprovided) data?
>
>
>>
>> Best wishes,
>>
>> Michael
>>
>> [[alternative HTML version deleted]]
>>
>> ______________________________________________
>> R-help at 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.
>
> David Winsemius
> Alameda, CA, USA
>
> [[alternative HTML version deleted]]
>
> ______________________________________________
> R-help at 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.
>
More information about the R-help
mailing list