|
|
A collection of R functions allowing probability simplexes to be drawn from collections of probability distributions.
CODE FOR SIMPLEX TRIANGLE
Simfunc<-function(a,b){
x<- -b/2 + 1-a
y<-b*sqrt(3)/2
return(c(x,y))
}
Sim<-function(x,y,z,aa){
if (min(x)<0){stop("Probabilities cannot be negative")}
if (min(y)<0){stop("Probabilities cannot be negative")}
a<-length(x)
vec<-rep(0,a)
for (i in 1:a){
vec[i]<-x[i]+y[i]
}
if (max(vec)>1){stop("Distributions sum to more than one")}
x1<-c(0,1)
y1<-c(0,0)
plot(x1,y1,type="l",xlim=c(-0.1,1.1),ylim=c(-0.1,sqrt(3/4)+0.1),xlab="", ylab="",axes=FALSE)
x2<-c(0,0.5)
y2<-c(0,sqrt(3/4))
par(new=T)
plot(x2,y2,type="l",xlim=c(-0.1,1.1),ylim=c(-0.1,sqrt(3/4)+0.1),xlab="", ylab="",axes=FALSE)
x3<-c(0.5,1)
y3<-c(sqrt(3/4),0)
par(new=T)
plot(x3,y3,type="l",xlim=c(-0.1,1.1),ylim=c(-0.1,sqrt(3/4)+0.1),xlab="", ylab="",axes=FALSE)
for (i in 1:a){
par(new=T)
plot(Simfunc(x[i],y[i])[1],Simfunc(x[i],y[i])[2],xlim=c(-0.1,1.1),ylim=c(-0.1,sqrt(3/4)+0.1),xlab="", ylab="",main="Simplex",axes=FALSE,pch=aa,col=z)
}
text(-0.05,-0.05,"(1,0,0)")
text(1.05,-0.05,"(0,0,1)")
text(0.5,0.92,"(0,1,0)")
}
Simline<-function(x,y,z){
if (min(x)<0){stop("Probabilities cannot be negative")}
if (min(y)<0){stop("Probabilities cannot be negative")}
a<-length(x)
vec<-rep(0,a)
for (i in 1:a){
vec[i]<-x[i]+y[i]
}
if (max(vec)>1){stop("Distributions sum to more than one")}
x1<-c(0,1)
y1<-c(0,0)
plot(x1,y1,type="l",xlim=c(-0.1,1.1),ylim=c(-0.1,sqrt(3/4)+0.1),xlab="", ylab="",axes=FALSE)
x2<-c(0,0.5)
y2<-c(0,sqrt(3/4))
par(new=T)
plot(x2,y2,type="l",xlim=c(-0.1,1.1),ylim=c(-0.1,sqrt(3/4)+0.1),xlab="", ylab="",axes=FALSE)
x3<-c(0.5,1)
y3<-c(sqrt(3/4),0)
par(new=T)
plot(x3,y3,type="l",xlim=c(-0.1,1.1),ylim=c(-0.1,sqrt(3/4)+0.1),xlab="", ylab="",axes=FALSE)
par(new=T)
a<-length(x)
plot(Simfunc(x,y)[1:a],Simfunc(x,y)[(a+1):(2*a)],type="l",xlim=c(-0.1,1.1),ylim=c(-0.1,sqrt(3/4)+0.1),xlab="", ylab="",main="Simplex",axes=FALSE,col=z)
text(-0.05,-0.05,"(1,0,0)")
text(1.05,-0.05,"(0,0,1)")
text(0.5,0.92,"(0,1,0)")
}
CODE FOR SIMPLEX TETRAHEDRON
#requires package "scatterplot3d"
#2d simplex
simplex<-function(x,y,z){
f1<-1-x-y/2
f2<-(sqrt(3)/2)*y
return(c(f1,f2))
}
#generation of tetrahedron
library(scatterplot3d)
x1<-c(0.0,0.5,1.0)
y1<-c(0,sqrt(3)/2,0)
z1<-c(0,0,0)
scatterplot3d(x1,y1,z1,xlim=c(-0.1,1),ylim=c(-0.1,1),zlim=c(-0.1,1),type="l",col.axis="blue",xlab="",ylab="",zlab="",grid=F)
x2<-c(0,0.5,1)
y2<-c(0,tan(pi/6)/2,0)
z2<-c(0,sqrt(2/3),0)
par(new=T)
scatterplot3d(x2,y2,z2,xlim=c(-0.1,1),ylim=c(-0.1,1),zlim=c(-0.1,1),type="l",col.axis="blue",xlab="",ylab="",zlab="",grid=F)
x3<-c(0.5,0.5)
y3<-c(tan(pi/6)/2,sqrt(3)/2)
z3<-c(sqrt(2/3),0)
par(new=T)
scatterplot3d(x3,y3,z3,xlim=c(-0.1,1),ylim=c(-0.1,1),zlim=c(-0.1,1),type="l",col.axis="blue",xlab="",ylab="",zlab="",grid=F)
x4<-c(0,1)
y4<-c(0,0)
z4<-c(0,0)
par(new=T)
scatterplot3d(x4,y4,z4,xlim=c(-0.1,1),ylim=c(-0.1,1),zlim=c(-0.1,1),type="l",col.axis="blue",xlab="",ylab="",zlab="",grid=F)
#convert length-four probability distribution into 3D point
ThreeD<-function(a,b,c,d){
if(d<1){
if(d>0){
a1<-simplex(a/(1-d),b/(1-d),c(1-d))[1]*(1-d)
a2<-simplex(a/(1-d),b/(1-d),c(1-d))[2]*(1-d)
x1<-(1-d)/2
y1<-(1-d)/(2*sqrt(3))
c1<-1/2-x1
c2<-1/(2*sqrt(3))-y1
e1<-simplex(a/(1-d),b/(1-d),c(1-d))[1]*(1-d)+c1
e2<-simplex(a/(1-d),b/(1-d),c(1-d))[2]*(1-d)+c2
e3<-sqrt(2/3)*d}
else{
e1<-simplex(a,b,c)[1]
e2<-simplex(a,b,c)[2]
e3<-0
}
}
else{
e1<-simplex(1/3,1/3,1/3)[1]
e2<-simplex(1/3,1/3,1/3)[2]
e3<-sqrt(2/3)
return(c(e1,e2,e3))
#plot (a,b,c,d)
par(new=T)
scatterplot3d(ThreeD(a,b,c,d)[1],ThreeD(a,b,c,d)[2],ThreeD(a,b,c,d)[3],xlim=c(-0.1,1),ylim=c(-0.1,1),zlim=c(-0.1,1),col.axis="blue",xlab="",ylab="",zlab="",grid=F)
The project summary page you can find here.