Thank you Bob.
I already had it worked.
Best,
M.
On Wed, Jan 6, 2016 at 1:19 PM, Bob Aronoff <bobaron...@gmail.com> wrote:
> I am working to understand the same issues with my datasets.
> Adapting Dr Zeileis' posting I have written a humble function.
> It creates a two cluster beta mix on a vector of data.
> It seems to be working well on my datasets.
> You are welcome to try on yours.
>
> regards,
> Bob
>
> _____
>
> bi.modal.beta<-function (vals){
>
> # vals is a vector data points in range (0,1). There should not be any
> 0,1, or NA in vals
>
> library(betareg)
> library(flexmix)
>
> d <- data.frame(y = vals)
>
> #create bimodal beta model
> m <- betamix(y ~ 1 | 1, data = d, k = 2)
>
> #extract beta parameters (mean and precision) with anti-link functions
> mu <- plogis(coef(m)[,1])
> phi <- exp(coef(m)[,2])
>
> #convert mean & precision to alpha & beta
> a <- mu * phi
> b <- (1 - mu) * phi
>
> #report parameters
> print("alpha:")
> print( a)
> print("beta:")
> print(b)
>
> #plot in order to inspect result
> ys <- seq(0, 1, by = 0.01)
> p <- prior(m$flexmix)
> # p is equivalent to the percentage of each distribution
>
> hist(d$y, breaks = 0:25/25, freq = FALSE,
> main = "Bimodal Betamix", xlab = "values")
> lines(ys, p[1] * dbeta(ys, shape1 = a[1], shape2 = b[1]) , lwd = 2,
> col="red")
> lines(ys, p[2] * dbeta(ys, shape1 = a[2], shape2 = b[2]) , lwd =
> 2,col="blue")
> }
>
> ____
>
> ______________________________________________
> 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.
>
[[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.
