# jh 2017.09.14: comments welcome

n=13

k=6 # faces (of say a die)

p0=(k:1)/sum(k:1)  # alter the probabilities as you wish
p0=rep(1/k,k)

final.max = k*n

s.min=1
s.max=k
prev.s.range=s.min:s.max

p.prev=p0

n.support.points = k*n

par(mfrow = c(3,4),mar=rep(0.0001,4))

y0 = -1.1*max(p0)
y1=   1.1*max(p0)
dy=max(p0)/k

for (i in 2:n){
	plot(c(s.min,s.max),c(1.15*y0,1.1*y1), 
	col="white", xaxt="n", yaxt="n")
	#abline(h=0)
	nn=length(prev.s.range)
	segments(prev.s.range,y0,
	         prev.s.range,y0+p.prev[prev.s.range],
	         col=prev.s.range+1,lwd=2)
	p=rep(0,final.max)
    txt = " "
	if(i==2) { txt = bquote(Y[1]) }
	if(i==3) { txt = bquote(Y[1]+Y[2]) }
	if(i==4) { txt = bquote(Y[1]+Y[2]+Y[3]) }
	if(i==5) { txt = bquote(Y[1]+....+Y[4]) }
	if(i==6) { txt = bquote(Y[1]+....+Y[5]) }
	if(i==7) { txt = bquote(Y[1]+....+Y[6]) }
	if(i==8) { txt = bquote(Y[1]+....+Y[7]) }
	if(i==9) { txt = bquote(Y[1]+....+Y[8]) }
	if(i==10) { txt = bquote(Y[1]+....+Y[9]) }
	if(i==11) { txt = bquote(Y[1]+....+Y[10]) }
	if(i==12) { txt = bquote(Y[1]+....+Y[11]) }
	if(i==13) { txt = bquote(Y[1]+....+Y[12]) }
	
	mass=p.prev[prev.s.range]
		m = sum(prev.s.range*mass)
		d = prev.s.range - m
		sd = sqrt(sum((d^2)*mass))
		skewness = sum(mass* ((prev.s.range-m)/sd)^3)
		kurtosis = sum(mass* ((prev.s.range-m)/sd)^4)
		Txt = toString(round(skewness,2))
		if(i<n) text(s.max,y0/2,Txt,adj=c(1,0.5))
		if(i==n) text(s.max,y0/2,
		    paste("Skewness:",Txt),
		  adj=c(1,0.5)) 
		Txt = toString(round(kurtosis,2))
		if(i<n) text(s.max,2*y0/3,Txt,
		   adj=c(1,0.5)) 
		if(i==n) text(s.max,2*y0/3,paste("Kurtosis:",Txt),
		   adj=c(1,0.5)) 



	text(mean(prev.s.range),y0,txt,adj=c(0.5,1.5),cex=1)
	for(prev.s in prev.s.range){
		if(prev.s == s.min | prev.s == s.max) text(prev.s,y0,
		 toString(prev.s),adj=c(0.5,1.5),cex=1)
		 
	    for (new in 1:k) {
			sum=prev.s+new
			prob = p.prev[prev.s] * p0[new]
			segments(prev.s,new*dy,
			         prev.s,new*dy+prob,
			         col=sum+1,lwd=2)
			p[sum] = p[sum] + prob
		}
	}
	s.min=i;
	s.max=i*k
    prev.s.range=s.min:s.max
	p.prev = p
}