### Test of the optimal Parameters routines.

set.seed(20160709)

## Set up model
pm1 <-
  TimelessNormalPM(muMean=c(Mechanics=2,Fluency=2),  varWeight=3,
                   precScale=solve(matrix(c(.7,.3,.3,.7),2,2)), Sdf=3)

## Draw parameters
p1 <- drawPMParam(pm1)

## Draw a fairly substantial sample.
stud1k <- drawInitialLatent(pm1,1000,param=p1)

## Used as starting point for iterative solutions.
p0 <- TimelessNormalParam(c(0,0),diag(2))


####MLE
p1.mle <- optimalPMParams(pm1,stud1k,Bayes=FALSE)
x2.mle <- sum((pvec(p1) - pvec(p1.mle$param))^2/abs(pvec(p1)))
stopifnot(x2.mle < qchisq(.99,length(pvec(p1))))

## Comparision to generic optimizer
go1.mle <- genericPMoptimizer(pm1,stud1k,Bayes=FALSE,
                              param=p0)
stopifnot (abs(pvec(p1.mle$param)-pvec(go1.mle$param))<.005,
           abs(p1.mle$dev -go1.mle$dev) <.005)


####MAP
p1.map <- optimalPMParams(pm1,stud1k,Bayes=TRUE)
x2.map <- sum((pvec(p1) - pvec(p1.map$param))^2/abs(pvec(p1)))
stopifnot(x2.map < qchisq(.99,length(pvec(p1))))
## Compare with generic optimizer.
go1.map <- genericPMoptimizer(pm1,stud1k,Bayes=TRUE,
                              param=p0)
stopifnot (abs(pvec(p1.map$param)-pvec(go1.map$param))<.005,
           abs(p1.map$dev -go1.map$dev) <.005)



p2a <- TimelessNormalParam(c(1,1),diag(2))
p2b <- TimelessNormalParam(c(0,0),diag(2))

## Draw a fairly substantial sample.
stud2a <- drawInitialLatent(pm1,500,param=p2a)
stud2b <- drawInitialLatent(pm1,500,param=p2b)
stud2 <- rbind(stud2a,stud2b)

## Equal weight, expect the mean around c(.5,.5)
p2.map <- optimalPMParams(pm1,stud2,weight=rep(.5,1000),Bayes=TRUE)
stopifnot(sum((getmu(p2.map$param)-c(.5,.5))^2)*500 < qchisq(.99,2))


go2.map <- genericPMoptimizer(pm1,stud2,weight=rep(.5,1000),Bayes=TRUE,
                              param=p0)
stopifnot (abs(pvec(p2.map$param)-pvec(go2.map$param))<.01,
           abs(p2.map$dev -go2.map$dev) <.01)


## 4 times the weight to the first data set, expect the mean around c(.8,.8)
p2w.map <- optimalPMParams(pm1,stud2,weights=rep(c(.8,.2),each=500),Bayes=TRUE)
stopifnot(sum((getmu(p2w.map)-c(.8,.8))^2)*500 < qchisq(.99,2))

go2w.map <- genericPMoptimizer(pm1,stud2,weights=rep(c(.8,.2),each=500),
                               Bayes=TRUE, param=p0)
stopifnot (abs(pvec(p2.map$param)-pvec(go2.map$param))<.01,
           abs(p2.map$dev -go2.map$dev) <.01)


p2a.map <- optimalPMParams(pm1,stud2a,Bayes=TRUE)$param
p2a2.map <- optimalPMParams(pm1,stud2a,2,Bayes=TRUE)$param
stopifnot(sum((getSigma(p2a.map)-getSigma(p2a))^2) < qchisq(.99,4))
stopifnot(sum((getSigma(p2a2.map)-getSigma(p2a))^2) < qchisq(.99,4))

go2a.map <- genericPMoptimizer(pm1,stud2a,Bayes=TRUE,param=p0)$param
go2a2.map <- genericPMoptimizer(pm1,stud2a,2,Bayes=TRUE,param=p0)$param
stopifnot(sum((getSigma(go2a.map)-getSigma(p2a))^2) < qchisq(.99,4))
stopifnot(sum((getSigma(go2a2.map)-getSigma(p2a))^2) < qchisq(.99,4))


p2a.mle <- optimalPMParams(pm1,stud2a,Bayes=FALSE)$param
p2a2.mle <- optimalPMParams(pm1,stud2a,2,Bayes=FALSE)$param
## These should be almost the same (minor difference as dividing by sum of
## the weights -1)
stopifnot(all(abs(getSigma(p2a.mle)-getSigma(p2a2.mle)) < .005))


go2a.mle <- optimalPMParams(pm1,stud2a,Bayes=FALSE)$param
go2a2.mle <- optimalPMParams(pm1,stud2a,2,Bayes=FALSE)$param
## These should be almost the same (minor difference as dividing by sum of
## the weights -1)
stopifnot(all(abs(getSigma(go2a.mle)-getSigma(go2a2.mle)) < .005))