Tag Archives: simplex

Confidence Regions for Parameters in the Simplex

Consider here the case where, in some parametric inference problem, parameter  is a point in the Simplex,

For instance, consider some regression, on compositional data,

> library(compositions)
>  data(DiagnosticProb)
>  Y=DiagnosticProb[,"type"]-1
>  X=DiagnosticProb[,c("A","B","C")]
>  model = glm(Y~ilr(X),family=binomial)
>  b = ilrInv(coef(model)[-1],orig=X)
>  as.numeric(b)
[1] 0.3447106 0.2374977 0.4177917

We can visualize that estimator on the simplex, using

>  tripoint=function(s){
+    p=s/sum(s)
+    abc2xy(matrix(p,1,3))
+  }

>  lab=LETTERS[1:3]
>  xl=c(-.1,1.25)
>  yl=c(-.1,1.15)
>  library(trifield)
>  A=abc2xy(matrix(c(1,0,0),1,3)) 
>  B=abc2xy(matrix(c(0,1,0),1,3))
>  C=abc2xy(matrix(c(0,0,1),1,3)) 
>  plot(0:1,0:1,col="white",
+  xlim=xl,ylim=yl,xlab="",ylab="",axes=FALSE)
>  polygon(rbind(A,B,C),col="light yellow")
>  text(B[1],-.05,lab[2])
>  text(A[1],1.05,lab[1])
>  text(C[1],-.05,lab[3])
>  segments((A[1]+C[1])/2,(A[2]+C[2])/2,B[1],B[2],col="grey",lty=2)
>  segments((A[1]+B[1])/2,(A[2]+B[2])/2,C[1],C[2],col="grey",lty=2)
>  segments((B[1]+C[1])/2,(B[2]+C[2])/2,A[1],A[2],col="grey",lty=2)
>  points(tripoint(b),pch=19,cex=2,col="red")

If we want to compute a ‘confidence region’, we can either use Bayesian models (with a Dirichlet distribution as prior distribution), or use bootstrap. We will use here the second idea

>  MB=matrix(NA,1e4,2)
>  for(sim in 1:1e4){
+    idx=sample(1:nrow(DiagnosticProb),
+    size=nrow(DiagnosticProb),replace=TRUE)
+  Y=DiagnosticProb[idx,"type"]-1
+  X=DiagnosticProb[idx,c("A","B","C")]
+  model = glm(Y~ilr(X),family=binomial)
+  MB[sim,]=tripoint(as.numeric(
+    ilrInv(coef(model)[-1],orig=X)))}

To get some ‘confidence region’, we can then use the bagplot, to get either a region where 50% of the boostraped estimators are, or 95%,

>  library(aplpack)
> P1=bagplot(MB[,1],MB[,2], factor =1.96, cex=.9,
+ dkmethod=2,show.baghull=TRUE) 
> P2=bagplot(MB[,1],MB[,2], factor =0.67, cex=.9,
+ dkmethod=2,show.baghull=TRUE) 

Then we can easily plot those two regions,

>  plot(0:1,0:1,col="white")
>  polygon(rbind(A,B,C),col="light yellow")
>  text(B[1],-.05,lab[2])
>  text(A[1],1.05,lab[1])
>  text(C[1],-.05,lab[3])
>  polygon(P1$hull.loop,col="yellow",border=NA)
>  polygon(P2$hull.loop,col="orange",border=NA)
>  segments((A[1]+C[1])/2,(A[2]+C[2])/2,B[1],B[2],col="grey",lty=2)
>  segments((A[1]+B[1])/2,(A[2]+B[2])/2,C[1],C[2],col="grey",lty=2)
>  segments((B[1]+C[1])/2,(B[2]+C[2])/2,A[1],A[2],col="grey",lty=2)
>  points(tripoint(b),pch=19,cex=2,col="red")

 

Overview on Multivariate Distributions

In June 2016, with Olivier L’Haridon, we will organize a (small) conference, in Rennes, on risk models in a multi-attribute framework. In order to fully enjoy the workshop (more to come on the blog), we will organize every month an internal workshop on that topic. We will start tomorrow afternoon, 13:00-14:30, and I will give a brief talk on multivariate distributions, with an emphasis on spherical / elliptical distributions, distributions on the simplex, and copulas. Slides are now online,

the Dirichlet distribution

In the course, since we are still introducing some concepts of dependent distributions, we will talk about the Dirichlet distribution, which is a distribution over the simplex of https://f.hypotheses.org/wp-content/blogs.dir/253/files/2017/07/diri11.gif. Let https://f.hypotheses.org/wp-content/blogs.dir/253/files/2017/07/diri01.gif denote the Gamma distribution with density (on https://f.hypotheses.org/wp-content/blogs.dir/253/files/2017/07/diri03.gif)

https://f.hypotheses.org/wp-content/blogs.dir/253/files/2017/07/diri02.gif

Let https://f.hypotheses.org/wp-content/blogs.dir/253/files/2017/07/diri04.gif denote independent https://f.hypotheses.org/wp-content/blogs.dir/253/files/2017/07/diri05.gif random variables, with https://f.hypotheses.org/wp-content/blogs.dir/253/files/2017/07/diri06.gif. Then https://f.hypotheses.org/wp-content/blogs.dir/253/files/2017/07/diri07.gif where

https://f.hypotheses.org/wp-content/blogs.dir/253/files/2017/07/diri08.gif

has a Dirichlet distribution with parameter

https://f.hypotheses.org/wp-content/blogs.dir/253/files/2017/07/diri09.gif

Note that https://f.hypotheses.org/wp-content/blogs.dir/253/files/2017/07/diri10.gif has a distribution in the simplex of https://f.hypotheses.org/wp-content/blogs.dir/253/files/2017/07/diri11.gif,

https://f.hypotheses.org/wp-content/blogs.dir/253/files/2017/07/diri40.gif

and has density

https://f.hypotheses.org/wp-content/blogs.dir/253/files/2017/07/diri12.gif

We will write https://f.hypotheses.org/wp-content/blogs.dir/253/files/2017/07/diri13.gif.

The density for different values of https://f.hypotheses.org/wp-content/blogs.dir/253/files/2017/07/diri20.gif can be visualized below, e.g. https://f.hypotheses.org/wp-content/blogs.dir/253/files/2017/07/diri21.gif, with some kind of symmetry,
https://f.hypotheses.org/wp-content/blogs.dir/253/files/2017/07/dirichlet222.gif
or https://f.hypotheses.org/wp-content/blogs.dir/253/files/2017/07/diri22.gif and https://f.hypotheses.org/wp-content/blogs.dir/253/files/2017/07/diri23.gif, below
https://f.hypotheses.org/wp-content/blogs.dir/253/files/2017/07/dirichlet522.gif
and finally, below, https://f.hypotheses.org/wp-content/blogs.dir/253/files/2017/07/diri24.gif


Note that marginal distributions are also Dirichlet, in the sense that if

https://f.hypotheses.org/wp-content/blogs.dir/253/files/2017/07/diri13.gif

then

https://f.hypotheses.org/wp-content/blogs.dir/253/files/2017/07/diri14.gif

if https://f.hypotheses.org/wp-content/blogs.dir/253/files/2017/07/diri15.gif, and if https://f.hypotheses.org/wp-content/blogs.dir/253/files/2017/07/diri16.gif, then https://f.hypotheses.org/wp-content/blogs.dir/253/files/2017/07/diri17.gif‘s have Beta distributions,

https://f.hypotheses.org/wp-content/blogs.dir/253/files/2017/07/diri18.gif

See Devroye (1986) section XI.4, or Frigyik, Kapila & Gupta (2010) .This distribution might also be called multivariate Beta distribution. In R, this function can be used as follows

> library(MCMCpack)
> alpha=c(2,2,5)
> x=seq(0,1,by=.05)
> vx=rep(x,length(x))
> vy=rep(x,each=length(x))
> vz=1-x-vy
> V=cbind(vx,vy,vz)
> D=ddirichlet(V, alpha)
> persp(x,x,matrix(D,length(x),length(x))

(to plot the density, as figures above). Note that we will come back on that distribution later on so-called Liouville copulas (see also Gupta & Richards (1986)).