# 12-15-97



go <- function() {
	mix.fig4 <- scan("fig4.data")
	xx <- seq(-4,4,l=201)
	ytrue <- .5*dnorm(xx,-1.5,sqrt(.5)) + .5*dnorm(xx,1,sqrt(1))
par(mfrow=c(2,2))
	fig4( mix.fig4, tk=seq(-4,4,length=3),head="m=2" )
		lines(xx,ytrue,lty=2)
	fig4( mix.fig4, tk=seq(-4,4,length=4),head="m=3" )
		lines(xx,ytrue,lty=2)
	fig4( mix.fig4, tk=seq(-4,4,length=5),head="m=4" )
		lines(xx,ytrue,lty=2)
	fig4( mix.fig4, tk=seq(-4,4,length=6),head="m=5" )
		lines(xx,ytrue,lty=2)

	fig4( mix.fig4, tk=seq(-4,4,length= 7),head="m= 6" )
		lines(xx,ytrue,lty=2)
	fig4( mix.fig4, tk=seq(-4,4,length= 8),head="m= 7" )
		lines(xx,ytrue,lty=2)
	fig4( mix.fig4, tk=seq(-4,4,length= 9),head="m= 8" )
		lines(xx,ytrue,lty=2)
	fig4( mix.fig4, tk=seq(-4,4,length=10),head="m= 9" )
		lines(xx,ytrue,lty=2)

	fig4( mix.fig4, tk=seq(-4,4,length=11),head="m=10" )
		lines(xx,ytrue,lty=2)
	fig4( mix.fig4, tk=seq(-4,4,length=12),head="m=11" )
		lines(xx,ytrue,lty=2)
	fig4( mix.fig4, tk=seq(-4,4,length=13),head="m=12" )
		lines(xx,ytrue,lty=2)
	fig4( mix.fig4, tk=seq(-4,4,length=14),head="m=13" )
		lines(xx,ytrue,lty=2)
}


# 9-18-97  [[12-15-97  replacing with new one ...see paper subdirectory]]

# NEW 4th Figure for Sagae/Scott paper

fig4 <- function( x, m=10, tk,head="" ) {	# version of cbh.c2.m1.s2
	if(missing(tk)){ab<-nicerange(x,0.0001)
		tk <- seq(ab[1],ab[2],length=m+1) }
		else{ m <- length(tk)-1 }
	n <- length(x)
	kbin <- integer(n)
	for(k in seq(m)) { kbin[ seq(n)[ x>=tk[k] & x<tk[k+1] ] ] <- k }
	i <- max(table(kbin))
browser()
	plot(0,0,xlim=c(-4,4),ylim=c(0,.3),type="n",
		xlab="x",ylab="",main=head)
#	plot(range(tk),c(-.4,1.6)*i/(n*(diff(range(tk))/m)),type="n")	# plot space

	mk0s <- rep(0,m); mk1s <- mk0s			# bin moments (padded by zero bin)
	hk <- diff(tk)					# bin widths
	mk <- 0.5 * (tk[-1]+tk[-(m+1)])			# bin centers

	for(k in seq(m)) {
	  xk <- x[kbin==k]			# data in bin
	  s0 <- length(xk); s1 <- 0
	  if(s0>0) { xk <- xk - mk[k] }		# centered on bin center
	  if(s0>0) { s0 <- s0/n; s1 <- sum(xk)/n }
	  mk0s[k] <- s0; mk1s[k] <- s1
	}

			# y <- c( c1...f1,c2...f2,...,cm...fm)
	A <- matrix(0,4*m,4*m)			# quadratic criterion
	B <- matrix(0,3*(m-1),4*m)		# continuity constraint matrix
	C <- rep(0,3*(m-1))			# right hand side constraints
	for(k in seq(1,length=m-1)) {		# set up equations
	  io <- 0				# row offset:  f continuous conditions
	  jo <- 4*(k-1)				# column offset: variables c,d,e,f
	  B[io+k,jo+1] <- hk[k]^2/12; B[io+k,jo+5] <- -hk[k+1]^2/12
	  B[io+k,jo+2] <- hk[k]^3/120; B[io+k,jo+6] <- hk[k+1]^3/120
	  B[io+k,jo+3] <- hk[k]^4/480; B[io+k,jo+7] <- -hk[k+1]^4/480
	  B[io+k,jo+4] <- hk[k]^5/6720; B[io+k,jo+8] <- hk[k+1]^5/6720
	  C[io+k] <- -mk0s[k]/hk[k] + mk0s[k+1]/hk[k+1] - 6*mk1s[k]/hk[k]^2 - 6*mk1s[k+1]/hk[k+1]^2
	  
	  io <- m-1				# row offset: f' continuous conditions
	  B[io+k,jo+1] <- hk[k]/2; B[io+k,jo+5] <- hk[k+1]/2
	  B[io+k,jo+2] <- hk[k]^2/10; B[io+k,jo+6] <- -hk[k+1]^2/10
	  B[io+k,jo+3] <- hk[k]^3/48; B[io+k,jo+7] <- hk[k+1]^3/48
	  B[io+k,jo+4] <- hk[k]^4/420; B[io+k,jo+8] <- -hk[k+1]^4/420
	  C[io+k] <- -12*mk1s[k]/hk[k]^3 + 12*mk1s[k+1]/hk[k+1]^3

	  io <- 2*(m-1)				# row offset: f'' continuous conditions
	  B[io+k,jo+1] <- 1; B[io+k,jo+5] <- -1
	  B[io+k,jo+2] <- hk[k]/2; B[io+k,jo+6] <- hk[k+1]/2
	  B[io+k,jo+3] <- hk[k]^2/8; B[io+k,jo+7] <- -hk[k+1]^2/8
	  B[io+k,jo+4] <- hk[k]^3/48; B[io+k,jo+8] <- hk[k+1]^3/48
	}

	for(k in seq(1,length=m)) {		# set up equations
	  jo <- 4*(k-1)				# column offset: variables c,d,e,f

	  A[jo+1,jo+1] <- hk[k]					# ck^2
	  A[jo+2,jo+2] <- hk[k]^3/12				# dk^2
	  A[jo+1,jo+3] <- A[jo+3,jo+1] <- hk[k]^3/12/2		# ck ek (split ck*ek, ek*ck)
	  A[jo+3,jo+3] <- hk[k]^5/320				# ek^2
	  A[jo+2,jo+4] <- A[jo+4,jo+2] <- hk[k]^5/240/2		# dk*fk  fk*dk
	  A[jo+4,jo+4] <- hk[k]^7/16128				# fk^2
	}
# browser()
	y <- solve(A) %*% t(B) %*% solve( B %*% solve(A) %*% t(B) ) %*% C
	i <- seq(0,by=4,length=m)
	ck <- y[1+i]; dk <- y[2+i]; ek <- y[3+i]; fk <- y[4+i]
	ak <- -ck*hk^2/24 - ek*hk^4/1920 + mk0s/hk
	bk <- -dk*hk^2/40 - fk*hk^4/4480 + 12*mk1s/hk^3

	for(k in seq(m)) {
	  tt <- seq(tk[k],tk[k+1],l=251)
	  tm <- tt - mk[k]			# centered
	  yy <- ak[k] + bk[k]*tm + ck[k]*tm^2/2 + dk[k]*tm^3/6 + ek[k]*tm^4/24 + fk[k]*tm^5/120
	  lines( tt , yy )
# browser()
	}
# browser()
# invisible(list(tk=tk,n=n,i=i,mk0s=mk0s,mk1s=mk1s,mks=mks,h=h))
}