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" <dwinsem...@comcast.net> wrote:


On Mar 21, 2017, at 5:04 AM, Michael Dayan <mdayan.resea...@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@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@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.


______________________________________________
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.

Reply via email to