Nothing
globalVariables('priors')
BiJGQD.mcmc=function(X,time,mesh=10,theta,sds,updates=10,burns=min(round(updates/2),25000),exclude=NULL,plot.chain=TRUE,RK.order=4,wrt=FALSE,Tag=NA,Dtype='Saddlepoint',Jdist='MVNormal',Jtype='Add',adapt=0,print.output=TRUE, decode=TRUE, palette = 'mono')
{
Mstar3 =1
recycle=FALSE
rtf=runif(2)
solver =function(Xs, Xt, theta, N , delt , N2, tt , P , alpha, lower , upper, tro ){}
rm(list =c('solver'))
theta = theta+runif(length(theta),0.001,0.002)*sign(theta)
T.seq=time
Dtypes =c('Saddlepoint','Normal')
Dindex = which(Dtypes==Dtype)
JDtypes=c('Normal','MVNormal')
JDindex = which(JDtypes==Jdist)
# A power function used to solve the syntax gap issue between R and C++
pow=function(x,p)
{
x^p
}
prod=function(a,b){a*b}
# Separate the time series into X and Y components
X1=X[,1]
X2=X[,2]
nnn=length(X1)
b1 = '\n==============================================================================\n'
b2 = '==============================================================================\n'
warn=c(
'1. Missing input: Argument {X} is missing.\n'
,'2. Missing input: Argument {time} is missing.\n'
,'3. Missing input: Argument {theta} is missing.\n'
,'4. Missing input: Argument {sds} is missing.\n'
,'5. Input type: Argument {X} must be of type matrix!.\n'
,'6. Input type: Argument {time} must be of type vector!.\n'
,'7. Input: Less starting parameters than model parameters.\n'
,'8. Input: More starting parameters than model parameters.\n'
,'9. Input: length(X) must be > 10.\n'
,'10. Input: length(time) must be > 10.\n'
,'11. Input: length(lower)!=1.\n'
,'12. Input: length(upper)!=1.\n'
,'13. Input: length(P)!=1.\n'
,'14. Input: length(mesh)!=1.\n'
,'15. Input: length(alpha)!=1.\n'
,'16. Input: length(Trunc)!=1.\n'
,'17. Input: length(RK.order)!=1.\n'
,'18. Density: Dtype has to be one of Saddlepoint or Normal.\n'
,'19. Density: Range [lower,upper] must be strictly positive for Dtype Gamma or InvGamma.\n'
,'20. Density: Dtype cannot be Beta for observations not in (0,1).\n'
,'21. Density: Argument {upper} must be > {lower}.\n'
,'22. Density: P must be >= 10.\n'
,'23. Density: Trunc[2] must be <= Trunc[1].\n'
,'24. ODEs : Large max(diff(time))/mesh may result in poor approximations. Try larger mesh.\n'
,'25. ODEs : max(diff(time))/mesh must be <1.\n'
,'26. ODEs : Runge Kutta scheme must be of order 4 or 10.\n'
,'27. ODEs : Argument {mesh} must be >= 5.\n'
,'28. Input: length(X)!=length(time).\n'
,'29. MCMC : Argument {burns} must be < {updates}.\n'
,'30. MCMC : Argument {updates} must be > 2.\n'
,'31. MCMC : length(theta)!=length(sds).\n'
,'32. Model: There has to be at least one model coefficient.\n'
,'33. Input: length(updates)!=1.\n'
,'34. Input: length(burns)!=1.\n'
,'35. Prior: priors(theta) must return a single value.\n'
,'36. Input: NAs not allowed.\n'
,'37. Input: length(Dtype)!=1.\n'
,'38. Input: NAs not allowed.\n'
,'39. Input: {Jdist} has to be of type "MVNormal".\n'
,'40. Input: {Jtype} has to be of type "Add" or "Mult".\n'
)
warntrue = rep(F,50)
check.thetas = function(theta,tt)
{
t=tt
theta = theta+runif(length(theta),0.001,0.002)*sign(theta)
#namess=c('a00','a10','a20','a01','a02','a11',
# 'b00','b10','b20','b01','b02','b11',
# 'c00','c10','c20','c01','c02','c11',
# 'd00','d10','d20','d01','d02','d11',
# 'e00','e10','e20','e01','e02','e11',
# 'f00','f10','f20','f01','f02','f11',
# 'Lam00','Lamy0','Nmu11','Nsig11','Nmu21','Nsig21','Nmu12','Nsig12','Nmu22','Nsig22','Lam10','Lamx2','Lamy1','Lamy2','Jmu1','Jmu2','Jsig11','Jsig12','Jsig22')
namess=c('a00','a10','a20','a01','a02','a11',
'b00','b10','b20','b01','b02','b11',
'c00','c10','c20','c01','c02','c11',
'd00','d10','d20','d01','d02','d11',
'e00','e10','e20','e01','e02','e11',
'f00','f10','f20','f01','f02','f11',
'Lam00','Lam10','Lam01','Jmu1','Jmu2','Jsig11','Jsig12','Jsig22')
func.list=rep(0,length(namess))
obs=objects(pos=1)
for(i in 1:length(namess))
{
if(sum(obs==namess[i])){func.list[i]=1}
}
pers.represented = rep(0,length(theta))
for(i in which(func.list==1))
{
for(j in 1:length(theta))
{
dresult1=eval(body(namess[i]))
theta[j] = theta[j]+runif(1,0.1,1)
dresult2=eval(body(namess[i]))
dff = abs(dresult1-dresult2)
if(any(round(dff,6)!=0)){pers.represented[j]=pers.represented[j]+1}
}
}
return(pers.represented)
}
check.thetas2 = function(theta)
{
namess=c('a00','a10','a20','a01','a02','a11',
'b00','b10','b20','b01','b02','b11',
'c00','c10','c20','c01','c02','c11',
'd00','d10','d20','d01','d02','d11',
'e00','e10','e20','e01','e02','e11',
'f00','f10','f20','f01','f02','f11',
'Lam00','Lam10','Lam01','Jmu1','Jmu2','Jsig11','Jsig12','Jsig22')
func.list=rep(0,length(namess))
obs=objects(pos=1)
for(i in 1:length(namess))
{
if(sum(obs==namess[i])){func.list[i]=1}
}
l=0
for(k in which(func.list==1))
{
str=body(namess[k])[2]
for(i in 1:length(theta))
{
for(j in 1:20)
{
str=sub(paste0('theta\\[',i,'\\]'),'clear',str)
}
}
l=l+length(grep('theta',str))
l
}
return(l)
}
# Check missing values first:
if(missing(X)) {warntrue[1]=T}
if(missing(time)) {warntrue[2]=T}
if(missing(theta)) {warntrue[3]=T}
if(missing(sds)) {warntrue[4]=T}
if(!is.matrix(X)) {warntrue[5]=T}
if(!is.vector(time)) {warntrue[6]=T}
# Check model parameters:
if(check.thetas2(theta)!=0) {warntrue[7]=T}
if(!warntrue[7]){if(any(check.thetas(theta,T.seq)==0)) {warntrue[8]=T}}
# Check input length:
if(dim(X)[1]<10) {warntrue[9]=T}
if(length(time)<10) {warntrue[10]=T}
if(length(mesh)!=1) {warntrue[14]=T}
if(length(RK.order)!=1) {warntrue[17]=T}
if(length(updates)!=1) {warntrue[33]=T}
if(length(burns)!=1) {warntrue[34]=T}
if(length(Dtype)!=1) {warntrue[37]=T}
# Check density approx parameters:
if(sum(Dindex)==0) {warntrue[18] =T}
# Miscelaneous checks:
excl=0
if(is.null(exclude)){excl=length(T.seq)-1+200}
if(!is.null(exclude)){excl=exclude}
test.this =max(diff(T.seq)[-excl])/mesh
if(test.this>0.1) {warntrue[24]=T}
if(test.this>=1) {warntrue[25]=T}
if(!warntrue[17]){if(!((RK.order==4)|(RK.order==10))) {warntrue[26]=T}}
if(!warntrue[14]){if(mesh<5) {warntrue[27]=T}}
if(dim(X)[1]!=length(time)) {warntrue[28]=T}
if(!any(warntrue[c(33,34)])){if(burns>updates) {warntrue[29]=T}}
if(!warntrue[33]){if(updates<2) {warntrue[30]=T}}
if(length(theta)!=length(sds)) {warntrue[31]=T}
if(any(is.na(X))||any(is.na(time))) {warntrue[36]=T}
if(sum(JDindex)==0) {warntrue[39] =TRUE}
if((Jtype!='Add')&&(Jtype!='Mult')) {warntrue[40] =TRUE}
# Print output:
if(any(warntrue))
{
prnt = b1
for(i in which(warntrue))
{
prnt = paste0(prnt,warn[i])
}
prnt = paste0(prnt,b2)
stop(prnt)
}
############################################################################
############################################################################
############################################################################
#==============================================================================
# Data resolution Module
#==============================================================================
#t=seq(0,100,by=1/100)
homo=T
homo.res=T
delt=(diff(T.seq)/mesh)[1]
t=T.seq
if(is.null(exclude))
{
if(sum(round(diff(T.seq)-diff(T.seq)[1],10)==0)!=length(T.seq)-1){homo.res=F}
}
if(!is.null(exclude))
{
if(sum(round(diff(T.seq)[-excl]-c(diff(T.seq)[-excl])[1],10)==0)!=length(T.seq)-1-length(excl)){homo.res=F}
}
#==============================================================================
# Prior distribution Module
#==============================================================================
if(sum(objects(pos=1)=='priors')==1)
{
pp=function(theta){}
body(pp)=parse(text =body(priors)[2])
prior.list=paste0('d(theta)',':',paste0(body(priors)[2]))
if(length(priors(theta))!=1){stop(" ==============================================================================
Incorrect input: Prior distribution must return a single value only!
==============================================================================");}
}
if(sum(objects(pos=1)=='priors')==0)
{
prior.list=paste0('d(theta)',':',' None.')
pp=function(theta){1}
}
if(length(pp(theta))!=1){stop("Prior density must return only a single value!")}
#==============================================================================
# Solution TYPE Module
#==============================================================================
namess=c('a00','a10','a20','a01','a02','a11',
'b00','b10','b20','b01','b02','b11',
'c00','c10','c20','c01','c02','c11',
'd00','d10','d20','d01','d02','d11',
'e00','e10','e20','e01','e02','e11',
'f00','f10','f20','f01','f02','f11')
func.list=rep(0,36)
obs=objects(pos=1)
for(i in 1:36)
{
if(sum(obs==namess[i])){func.list[i]=1}
}
set1.1=c(3,14,15,20,21) # X Nonlinear
set1.2=c(11,28,29,34,35) # Y Nonlinear
set4=c(3,5,6,9,11,12,14:18,20:24,26:30,32:36) # Nonlinear terms
set5=c(4:6,8,9,12,16:18,22:24,26,27,30,32,33,36) # Dependant terms
state3.0=( any(func.list[set1.1]==1) & any(func.list[set1.2]==1) &(sum(func.list[set5])==0)) # Any nonlinear in BOTH and NO interaction term
state3.1=((sum(any(func.list[set1.1]==1))==0)&( any(func.list[set1.2]==1) )&(sum(func.list[set5])==0)) # Any nonlinear Y ,Normal X and NO interaction term
state3.2=( any(func.list[set1.1]==1) &(sum(any(func.list[set1.2]==1))==0)&(sum(func.list[set5])==0)) # Any nonlinear X ,Normal Y and NO interaction term
state4 =(any(func.list[set4]==1)&any(func.list[set5]==1)) # Any nonlinear term AND interaction terms
state1 =!(state3.0|state3.1|state3.2|state4)
state1=F
state3.0=F
state3.1=F
state3.2=F
state4=T
sol.state=c(state1,state3.0,state3.1,state3.2,state4)
sol.state=c('2nd Ord. Truncation, Bivariate Normal','4th Ord. Truncation, Independent Saddlepoints ','2nd & 4th Ord. Truncation (Indep)','4th & 2nd Ord. Trunc. (Indep)',"4th Ord. Truncation, Bivariate-Saddlepoint")[which(sol.state==T)]
#==============================================================================
# Check the syntax of each coefficient
#==============================================================================
strcheck=function(str)
{
k=0
if(length(grep('/',str))!=0)
{
warning('C++ uses integer division when denominators are of type int. Use decimal values if possible (e.g. 0.5 i.s.o. 1/2).',call. = F)
k=k+1
}
if(length(grep('\\^',str))!=0)
{
stop('^-Operator not overloaded in C++: Use pow(x,p) instead (e.g. pow(x,2) i.s.o. x^2).',call. = F)
k=k+1
}
return(k)
}
for(i in which(func.list==1))
{
strcheck(body(namess[i])[2])
}
############################################################################
############################################################################
############################################################################
tro=c(5,8,8,8,14)[which(c(state1,state3.0,state3.1,state3.2,state4)==T)]
tro = 14
if(tro==14)
{
txtA='
#include <RcppArmadillo.h>
#include <math.h>
#include <Rcpp.h>
#define pi 3.14159265358979323846 /* pi */
using namespace arma;
using namespace Rcpp;
using namespace R;
// [[Rcpp::depends("RcppArmadillo")]]
// [[Rcpp::export]]
vec prod(vec a,vec b)
{
return(a%b);
}
mat f(mat a,vec theta,vec t,int N2)
{
mat atemp(N2,34);'
secmom = 5
}
if(RK.order==4)
{
if(homo.res)
{
txtB= '
return atemp;
}
// [[Rcpp::export]]
List solver(vec Xs,vec Ys,vec Xt,vec Yt,vec theta,int N,double delt,int N2,vec tt,mat starts,int tro,int secmom,vec seq1,vec seq2)
{
mat resss(N2,3);
mat x0(N2,tro);
mat xa(N2,tro);
mat xe(N2,tro);
mat fx1(N2,tro);
mat fx2(N2,tro);
mat fx3(N2,tro);
mat fx4(N2,tro);
mat fx5(N2,tro);
mat fx6(N2,tro);
double whch =0;
x0.fill(0);
for (int i = 1; i <= 14; i++)
{
x0.col(i-1)=pow(Xs,seq1[i-1])%pow(Ys,seq2[i-1]);
x0.col(i-1+14)=pow(Xs,seq1[i-1])%pow(Ys,seq2[i-1]);
}
vec d=tt;
for (int i = 1; i < N; i++)
{
fx1 = f(x0,theta,d,N2)*delt;
fx2 = f(x0+0.25*fx1,theta,d+0.25*delt,N2)*delt;
fx3 = f(x0+0.09375*fx1+0.28125*fx2,theta,d+0.375*delt,N2)*delt;
fx4 = f(x0+0.879381*fx1-3.277196*fx2+ 3.320892*fx3,theta,d+0.9230769*delt,N2)*delt;
fx5 = f(x0+2.032407*fx1-8*fx2+7.173489*fx3-0.2058967*fx4,theta,d+delt,N2)*delt;
fx6 = f(x0-0.2962963*fx1+2*fx2-1.381676*fx3+0.4529727*fx4-0.275*fx5,theta,d+0.5*delt,N2)*delt;
xa = x0+0.1185185*fx1+0.5189864*fx3+0.5061315*fx4-0.18*fx5+0.03636364*fx6;
x0 = x0+0.1157407*fx1+0.5489279*fx3+0.5353314*fx4-0.2*fx5;
xe = abs(x0.col(1)-xa.col(1));
if(xe.max()>whch)
{
whch = xe.max();
}
d=d+delt;
}
'
}
if(!homo.res)
{
txtB= '
return atemp;
}
// [[Rcpp::export]]
List solver(vec Xs,vec Ys,vec Xt,vec Yt,vec theta,int N, mat delt,int N2,vec d,mat starts,int tro,int secmom,vec seq1,vec seq2)
{
mat resss(N2,3);
mat x0(N2,tro);
mat xa(N2,tro);
mat xe(N2,tro);
mat fx1(N2,tro);
mat fx2(N2,tro);
mat fx3(N2,tro);
mat fx4(N2,tro);
mat fx5(N2,tro);
mat fx6(N2,tro);
double whch =0;
x0.fill(0);
for (int i = 1; i <= 14; i++)
{
x0.col(i-1)=pow(Xs,seq1[i-1])%pow(Ys,seq2[i-1]);
x0.col(i-1+14)=pow(Xs,seq1[i-1])%pow(Ys,seq2[i-1]);
}
for (int i = 1; i < N; i++)
{
fx1 = f(x0,theta,d,N2)%delt;
fx2 = f(x0+0.25*fx1,theta,d+0.25*delt.col(0),N2)%delt;
fx3 = f(x0+0.09375*fx1+0.28125*fx2,theta,d+0.375*delt.col(0),N2)%delt;
fx4 = f(x0+0.879381*fx1-3.277196*fx2+ 3.320892*fx3,theta,d+0.9230769*delt.col(0),N2)%delt;
fx5 = f(x0+2.032407*fx1-8*fx2+7.173489*fx3-0.2058967*fx4,theta,d+delt.col(0),N2)%delt;
fx6 = f(x0-0.2962963*fx1+2*fx2-1.381676*fx3+0.4529727*fx4-0.275*fx5,theta,d+0.5*delt.col(0),N2)%delt;
xa = x0+0.1185185*fx1+0.5189864*fx3+0.5061315*fx4-0.18*fx5+0.03636364*fx6;
x0 = x0+0.1157407*fx1+0.5489279*fx3+0.5353314*fx4-0.2*fx5;
xe = abs(x0.col(1)-xa.col(1));
if(xe.max()>whch)
{
whch = xe.max();
}
d=d+delt.col(0);
}
'
}
}
if(RK.order==10)
{
if(homo.res)
{
txtB= '
return atemp;
}
// [[Rcpp::export]]
List solver(vec Xs,vec Ys,vec Xt,vec Yt,vec theta,int N,double delt,int N2,vec tt,mat starts,int tro,int secmom,vec seq1,vec seq2)
{
mat resss(N2,3);
mat fx0(N2,tro);
mat fx1(N2,tro);
mat fx2(N2,tro);
mat fx3(N2,tro);
mat fx4(N2,tro);
mat fx5(N2,tro);
mat fx6(N2,tro);
mat fx7(N2,tro);
mat fx8(N2,tro);
mat fx9(N2,tro);
mat fx10(N2,tro);
mat fx11(N2,tro);
mat fx12(N2,tro);
mat fx13(N2,tro);
mat fx14(N2,tro);
mat fx15(N2,tro);
mat fx16(N2,tro);
mat x0(N2,tro);
mat x1(N2,tro);
mat x2(N2,tro);
mat x3(N2,tro);
mat x4(N2,tro);
mat x5(N2,tro);
mat x6(N2,tro);
mat x7(N2,tro);
mat x8(N2,tro);
mat x9(N2,tro);
mat x10(N2,tro);
mat x11(N2,tro);
mat x12(N2,tro);
mat x13(N2,tro);
mat x14(N2,tro);
mat x15(N2,tro);
mat x16(N2,tro);
x0.fill(0);
double whch =0;
for (int i = 1; i <= 14; i++)
{
x0.col(i-1)=pow(Xs,seq1[i-1])%pow(Ys,seq2[i-1]);
x0.col(i-1+14)=pow(Xs,seq1[i-1])%pow(Ys,seq2[i-1]);
}
vec d=tt;
for (int i = 1; i < N; i++)
{
fx0=f(x0,theta,d,N2);
x1=x0+delt*(0.1*fx0);
fx1=f(x1,theta,d+0.100000000000000000000000000000000000000000000000000000000000*delt,N2);
x2=x0+delt*(-0.915176561375291*fx0+1.45453440217827*fx1);
fx2=f(x2,theta,d+0.539357840802981787532485197881302436857273449701009015505500*delt,N2);
x3=x0+delt*(0.202259190301118*fx0+0.606777570903354*fx2);
fx3=f(x3,theta,d+0.809036761204472681298727796821953655285910174551513523258250*delt,N2);
x4=x0+delt*(0.184024714708644*fx0+0.197966831227192*fx2-0.0729547847313633*fx3);
fx4=f(x4,theta,d+0.309036761204472681298727796821953655285910174551513523258250*delt,N2);
x5=x0+delt*(0.0879007340206681*fx0+0.410459702520261*fx3+0.482713753678866*fx4);
fx5=f(x5,theta,d+0.981074190219795268254879548310562080489056746118724882027805*delt,N2);
x6=x0+delt*(0.085970050490246*fx0+0.330885963040722*fx3+0.48966295730945*fx4-0.0731856375070851*fx5);
fx6=f(x6,theta,d+0.833333333333333333333333333333333333333333333333333333333333*delt,N2);
x7=x0+delt*(0.120930449125334*fx0+0.260124675758296*fx4+0.0325402621549091*fx5-0.0595780211817361*fx6);
fx7=f(x7,theta,d+0.354017365856802376329264185948796742115824053807373968324184*delt,N2);
x8=x0+delt*(0.110854379580391*fx0-0.0605761488255006*fx5+0.321763705601778*fx6+0.510485725608063*fx7);
fx8=f(x8,theta,d+0.882527661964732346425501486979669075182867844268052119663791*delt,N2);
x9=x0+delt*(0.112054414752879*fx0-0.144942775902866*fx5-0.333269719096257*fx6+0.49926922955688*fx7+0.509504608929686*fx8);
fx9=f(x9,theta,d+0.642615758240322548157075497020439535959501736363212695909875*delt,N2);
x10=x0+delt*(0.113976783964186*fx0-0.0768813364203357*fx5+0.239527360324391*fx6+0.397774662368095*fx7+0.0107558956873607*fx8-0.327769124164019*fx9);
fx10=f(x10,theta,d+0.357384241759677451842924502979560464040498263636787304090125*delt,N2);
x11=x0+delt*(0.0798314528280196*fx0-0.0520329686800603*fx5-0.0576954146168549*fx6+0.194781915712104*fx7+0.145384923188325*fx8-0.0782942710351671*fx9-0.114503299361099*fx10);
fx11=f(x11,theta,d+0.117472338035267653574498513020330924817132155731947880336209*delt,N2);
x12=x0+delt*(0.985115610164857*fx0+0.330885963040722*fx3+0.48966295730945*fx4-1.37896486574844*fx5-0.861164195027636*fx6+5.78428813637537*fx7+3.28807761985104*fx8-2.38633905093136*fx9-3.25479342483644*fx10-2.16343541686423*fx11);
fx12=f(x12,theta,d+0.833333333333333333333333333333333333333333333333333333333333*delt,N2);
x13=x0+delt*(0.895080295771633*fx0+0.197966831227192*fx2-0.0729547847313633*fx3-0.851236239662008*fx5+0.398320112318533*fx6+3.63937263181036*fx7+1.5482287703983*fx8-2.12221714704054*fx9-1.58350398545326*fx10-1.71561608285936*fx11-0.0244036405750127*fx12);
fx13=f(x13,theta,d+0.309036761204472681298727796821953655285910174551513523258250*delt,N2);
x14=x0+delt*(-0.915176561375291*fx0+1.45453440217827*fx1+0*fx2+0*fx3-0.777333643644968*fx4+0*fx5-0.0910895662155176*fx6+0.0910895662155176*fx12+0.777333643644968*fx13);
fx14=f(x14,theta,d+0.539357840802981787532485197881302436857273449701009015505500*delt,N2);
x15=x0+delt*(0.1*fx0-0.157178665799771*fx2+0.157178665799771*fx14);
fx15=f(x15,theta,d+0.100000000000000000000000000000000000000000000000000000000000*delt,N2);
x16=x0+delt*(0.181781300700095*fx0+0.675*fx1+0.34275815984719*fx2+0*fx3+0.259111214548323*fx4-0.358278966717952*fx5-1.04594895940883*fx6+0.930327845415627*fx7+1.77950959431708*fx8+0.1*fx9-0.282547569539044*fx10-0.159327350119973*fx11-0.145515894647002*fx12-0.259111214548323*fx13-0.34275815984719*fx14-0.675*fx15);
fx16=f(x16,theta,d+delt,N2);
x0=x0+(0.0333333333333333333333333333333333333333333333333333333333333*fx0
+0.0250000000000000000000000000000000000000000000000000000000000*fx1
+0.0333333333333333333333333333333333333333333333333333333333333*fx2
+0.000000000000000000000000000000000000000000000000000000000000*fx3
+0.0500000000000000000000000000000000000000000000000000000000000*fx4
+0.000000000000000000000000000000000000000000000000000000000000*fx5
+0.0400000000000000000000000000000000000000000000000000000000000*fx6
+0.000000000000000000000000000000000000000000000000000000000000*fx7
+0.189237478148923490158306404106012326238162346948625830327194*fx8
+0.277429188517743176508360262560654340428504319718040836339472*fx9
+0.277429188517743176508360262560654340428504319718040836339472*fx10
+0.189237478148923490158306404106012326238162346948625830327194*fx11
-0.0400000000000000000000000000000000000000000000000000000000000*fx12
-0.0500000000000000000000000000000000000000000000000000000000000*fx13
-0.0333333333333333333333333333333333333333333333333333333333333*fx14
-0.0250000000000000000000000000000000000000000000000000000000000*fx15
+0.0333333333333333333333333333333333333333333333333333333333333*fx16)*delt;
d=d+delt;
}
'
}
if(!homo.res)
{
txtB= '
return atemp;
}
// [[Rcpp::export]]
List solver(vec Xs,vec Ys,vec Xt,vec Yt,vec theta,int N,double delt,int N2,vec tt,mat starts,int tro,int secmom,vec seq1,vec seq2)
{
mat resss(N2,3);
mat fx0(N2,tro);
mat fx1(N2,tro);
mat fx2(N2,tro);
mat fx3(N2,tro);
mat fx4(N2,tro);
mat fx5(N2,tro);
mat fx6(N2,tro);
mat fx7(N2,tro);
mat fx8(N2,tro);
mat fx9(N2,tro);
mat fx10(N2,tro);
mat fx11(N2,tro);
mat fx12(N2,tro);
mat fx13(N2,tro);
mat fx14(N2,tro);
mat fx15(N2,tro);
mat fx16(N2,tro);
mat x0(N2,tro);
mat x1(N2,tro);
mat x2(N2,tro);
mat x3(N2,tro);
mat x4(N2,tro);
mat x5(N2,tro);
mat x6(N2,tro);
mat x7(N2,tro);
mat x8(N2,tro);
mat x9(N2,tro);
mat x10(N2,tro);
mat x11(N2,tro);
mat x12(N2,tro);
mat x13(N2,tro);
mat x14(N2,tro);
mat x15(N2,tro);
mat x16(N2,tro);
double whch =0;
x0.fill(0);
for (int i = 1; i <= 14; i++)
{
x0.col(i-1)=pow(Xs,seq1[i-1])%pow(Ys,seq2[i-1]);
x0.col(i-1+14)=pow(Xs,seq1[i-1])%pow(Ys,seq2[i-1]);
}
for (int i = 1; i < N; i++)
{
fx0=f(x0,theta,d,N2);
x1=x0+delt%(0.1*fx0);
fx1=f(x1,theta,d+0.100000000000000000000000000000000000000000000000000000000000*delt.col(0),N2);
x2=x0+delt%(-0.915176561375291*fx0+1.45453440217827*fx1);
fx2=f(x2,theta,d+0.539357840802981787532485197881302436857273449701009015505500*delt.col(0),N2);
x3=x0+delt%(0.202259190301118*fx0+0.606777570903354*fx2);
fx3=f(x3,theta,d+0.809036761204472681298727796821953655285910174551513523258250*delt.col(0),N2);
x4=x0+delt%(0.184024714708644*fx0+0.197966831227192*fx2-0.0729547847313633*fx3);
fx4=f(x4,theta,d+0.309036761204472681298727796821953655285910174551513523258250*delt.col(0),N2);
x5=x0+delt%(0.0879007340206681*fx0+0.410459702520261*fx3+0.482713753678866*fx4);
fx5=f(x5,theta,d+0.981074190219795268254879548310562080489056746118724882027805*delt.col(0),N2);
x6=x0+delt%(0.085970050490246*fx0+0.330885963040722*fx3+0.48966295730945*fx4-0.0731856375070851*fx5);
fx6=f(x6,theta,d+0.833333333333333333333333333333333333333333333333333333333333*delt.col(0),N2);
x7=x0+delt%(0.120930449125334*fx0+0.260124675758296*fx4+0.0325402621549091*fx5-0.0595780211817361*fx6);
fx7=f(x7,theta,d+0.354017365856802376329264185948796742115824053807373968324184*delt.col(0),N2);
x8=x0+delt%(0.110854379580391*fx0-0.0605761488255006*fx5+0.321763705601778*fx6+0.510485725608063*fx7);
fx8=f(x8,theta,d+0.882527661964732346425501486979669075182867844268052119663791*delt.col(0),N2);
x9=x0+delt%(0.112054414752879*fx0-0.144942775902866*fx5-0.333269719096257*fx6+0.49926922955688*fx7+0.509504608929686*fx8);
fx9=f(x9,theta,d+0.642615758240322548157075497020439535959501736363212695909875*delt.col(0),N2);
x10=x0+delt%(0.113976783964186*fx0-0.0768813364203357*fx5+0.239527360324391*fx6+0.397774662368095*fx7+0.0107558956873607*fx8-0.327769124164019*fx9);
fx10=f(x10,theta,d+0.357384241759677451842924502979560464040498263636787304090125*delt.col(0),N2);
x11=x0+delt%(0.0798314528280196*fx0-0.0520329686800603*fx5-0.0576954146168549*fx6+0.194781915712104*fx7+0.145384923188325*fx8-0.0782942710351671*fx9-0.114503299361099*fx10);
fx11=f(x11,theta,d+0.117472338035267653574498513020330924817132155731947880336209*delt.col(0),N2);
x12=x0+delt%(0.985115610164857*fx0+0.330885963040722*fx3+0.48966295730945*fx4-1.37896486574844*fx5-0.861164195027636*fx6+5.78428813637537*fx7+3.28807761985104*fx8-2.38633905093136*fx9-3.25479342483644*fx10-2.16343541686423*fx11);
fx12=f(x12,theta,d+0.833333333333333333333333333333333333333333333333333333333333*delt.col(0),N2);
x13=x0+delt%(0.895080295771633*fx0+0.197966831227192*fx2-0.0729547847313633*fx3-0.851236239662008*fx5+0.398320112318533*fx6+3.63937263181036*fx7+1.5482287703983*fx8-2.12221714704054*fx9-1.58350398545326*fx10-1.71561608285936*fx11-0.0244036405750127*fx12);
fx13=f(x13,theta,d+0.309036761204472681298727796821953655285910174551513523258250*delt.col(0),N2);
x14=x0+delt%(-0.915176561375291*fx0+1.45453440217827*fx1+0*fx2+0*fx3-0.777333643644968*fx4+0*fx5-0.0910895662155176*fx6+0.0910895662155176*fx12+0.777333643644968*fx13);
fx14=f(x14,theta,d+0.539357840802981787532485197881302436857273449701009015505500*delt.col(0),N2);
x15=x0+delt%(0.1*fx0-0.157178665799771*fx2+0.157178665799771*fx14);
fx15=f(x15,theta,d+0.100000000000000000000000000000000000000000000000000000000000*delt.col(0),N2);
x16=x0+delt%(0.181781300700095*fx0+0.675*fx1+0.34275815984719*fx2+0*fx3+0.259111214548323*fx4-0.358278966717952*fx5-1.04594895940883*fx6+0.930327845415627*fx7+1.77950959431708*fx8+0.1*fx9-0.282547569539044*fx10-0.159327350119973*fx11-0.145515894647002*fx12-0.259111214548323*fx13-0.34275815984719*fx14-0.675*fx15);
fx16=f(x16,theta,d+delt.col(0),N2);
x0=x0+(0.0333333333333333333333333333333333333333333333333333333333333*fx0
+0.0250000000000000000000000000000000000000000000000000000000000*fx1
+0.0333333333333333333333333333333333333333333333333333333333333*fx2
+0.000000000000000000000000000000000000000000000000000000000000*fx3
+0.0500000000000000000000000000000000000000000000000000000000000*fx4
+0.000000000000000000000000000000000000000000000000000000000000*fx5
+0.0400000000000000000000000000000000000000000000000000000000000*fx6
+0.000000000000000000000000000000000000000000000000000000000000*fx7
+0.189237478148923490158306404106012326238162346948625830327194*fx8
+0.277429188517743176508360262560654340428504319718040836339472*fx9
+0.277429188517743176508360262560654340428504319718040836339472*fx10
+0.189237478148923490158306404106012326238162346948625830327194*fx11
-0.0400000000000000000000000000000000000000000000000000000000000*fx12
-0.0500000000000000000000000000000000000000000000000000000000000*fx13
-0.0333333333333333333333333333333333333333333333333333333333333*fx14
-0.0250000000000000000000000000000000000000000000000000000000000*fx15
+0.0333333333333333333333333333333333333333333333333333333333333*fx16)%delt;
d=d+delt.col(0);
}
'
}
}
if(state4)
{
if(Dtype=='Saddlepoint')
{
txtC='
vec probs=exp(-x0.col(28)-x0.col(29)-x0.col(30)-x0.col(31)-x0.col(32)-x0.col(33));
vec m00 = (1+0*x0.col(0));
vec m10 = (x0.col(0) -probs%x0.col(14))/(1-probs);
vec m20 = (x0.col(1) -probs%x0.col(15))/(1-probs);
vec m30 = (x0.col(2) -probs%x0.col(16))/(1-probs);
vec m40 = (x0.col(3) -probs%x0.col(17))/(1-probs);
vec m01 = (x0.col(4) -probs%x0.col(18))/(1-probs);
vec m02 = (x0.col(5) -probs%x0.col(19))/(1-probs);
vec m03 = (x0.col(6) -probs%x0.col(20))/(1-probs);
vec m04 = (x0.col(7) -probs%x0.col(21))/(1-probs);
vec m11 = (x0.col(8) -probs%x0.col(22))/(1-probs);
vec m12 = (x0.col(9) -probs%x0.col(23))/(1-probs);
vec m21 = (x0.col(10)-probs%x0.col(24))/(1-probs);
vec m22 = (x0.col(11)-probs%x0.col(25))/(1-probs);
vec m13 = (x0.col(12)-probs%x0.col(26))/(1-probs);
vec m31 = (x0.col(13)-probs%x0.col(27))/(1-probs);
vec mm00 = (1+0*x0.col(0));
vec mm10 = x0.col(14);
vec mm20 = x0.col(15);
vec mm30 = x0.col(16);
vec mm40 = x0.col(17);
vec mm01 = x0.col(18);
vec mm02 = x0.col(19);
vec mm03 = x0.col(20);
vec mm04 = x0.col(21);
vec mm11 = x0.col(22);
vec mm12 = x0.col(23);
vec mm21 = x0.col(24);
vec mm22 = x0.col(25);
vec mm13 = x0.col(26);
vec mm31 = x0.col(27);
vec k10 = m10;
vec k20 = m20-pow(m10,2);
vec k30 = m30-3*m20%m10+2*pow(m10,3);
vec k40 = m40 -4*m30%m10-3*pow(m20,2)+12*m20%pow(m10,2)-6*pow(m10,4);
vec k01 = m01;
vec k02 = m02- pow(m01,2) ;
vec k03 = m03-3*m02%m01+2*pow(m01,3) ;
vec k04 = m04-4*m03%m01-3*pow(m02,2)+12*m02%pow(m01,2)-6*pow(m01,4);
vec k11 = m11-m10%m01;
vec k21 = m21-2*m11%m10-m20%m01+2*pow(m10,2)%m01;
vec k12 = m12-2*m11%m01-m02%m10+2*pow(m01,2)%m10;
vec k22 = m22-2*m21%m01-2*m12%m10-m20%m02-2*pow(m11,2)+8*m11%m01%m10+2*m02%pow(m10,2)+2*m20%pow(m01,2)-6*pow(m10,2)%pow(m01,2) ;
vec k31 = m31-3*m21%m10-m30%m01-3*m20%m11+6*m11%pow(m10,2)+6*m20%m10%m01-6*pow(m10,3)%m01 ;
vec k13 = m13-3*m12%m01-m03%m10-3*m02%m11+6*m11%pow(m01,2)+6*m02%m01%m10-6*pow(m01,3)%m10;
vec kk10 = mm10;
vec kk20 = mm20-pow(mm10,2);
vec kk30 = mm30-3*mm20%mm10+2*pow(mm10,3);
vec kk40 = mm40 -4*mm30%mm10-3*pow(mm20,2)+12*mm20%pow(mm10,2)-6*pow(mm10,4);
vec kk01 = mm01;
vec kk02 = mm02- pow(mm01,2) ;
vec kk03 = mm03-3*mm02%mm01+2*pow(mm01,3) ;
vec kk04 = mm04-4*mm03%mm01-3*pow(mm02,2)+12*mm02%pow(mm01,2)-6*pow(mm01,4);
vec kk11 = mm11-mm10%mm01;
vec kk21 = mm21-2*mm11%mm10-mm20%mm01+2*pow(mm10,2)%mm01;
vec kk12 = mm12-2*mm11%mm01-mm02%mm10+2*pow(mm01,2)%mm10;
vec kk22 = mm22-2*mm21%mm01-2*mm12%mm10-mm20%mm02-2*pow(mm11,2)+8*mm11%mm01%mm10+2*mm02%pow(mm10,2)+2*mm20%pow(mm01,2)-6*pow(mm10,2)%pow(mm01,2) ;
vec kk31 = mm31-3*mm21%mm10-mm30%mm01-3*mm20%mm11+6*mm11%pow(mm10,2)+6*mm20%mm10%mm01-6*pow(mm10,3)%mm01 ;
vec kk13 = mm13-3*mm12%mm01-mm03%mm10-3*mm02%mm11+6*mm11%pow(mm01,2)+6*mm02%mm01%mm10-6*pow(mm01,3)%mm10;
vec a(N2);
vec b(N2);
vec abser(N2);
abser=0.1+abser;
a.ones();
b.ones();
vec det=(k10%k01-k11%k11);
a=-(Xt-k10)%k20/det+(Yt-k01)%k11/det;
b=+(Xt-k10)%k11/det-(Yt-k01)%k02/det;
vec gg(N2);
vec hh(N2);
vec gg1(N2);
vec hh1(N2);
vec gg2(N2);
vec hh2(N2);
vec ar(N2);
vec br(N2);
vec anew(N2);
vec bnew(N2);
int ind=0;
while((max(abser)>0.001)&&(ind<1500))
{
gg=k10+k20%a+(1.0/2.0)*k30%a%a+(1.0/6.0)*k40%a%a%a +k11%b +(1.0/2.0)*k12%b%b+k21%a%b+(1.0/6.0)*b%b%b%k13+(1.0/2.0)*a%a%b%k31+(1.0/2.0)*a%b%b%k22-Xt;
hh=k01+k02%b+(1.0/2.0)*k03%b%b+(1.0/6.0)*k04%b%b%b +k11%a +k12%a%b+(1.0/2.0)*k21%a%a+(1.0/2.0)*a%b%b%k13+(1.0/6.0)*a%a%a%k31+(1.0/2.0)*a%a%b%k22-Yt;
gg1=k20+k30%a+(1.0/2.0)*k40%a%a+k21%b+a%b%k31+(1.0/2.0)*b%b%k22;
gg2=k11 +k12%b+k21%a+(1.0/2.0)*b%b%k13+(1.0/2.0)*a%a%k31+a%b%k22;
hh1=k11 +k12%b+k21%a+(1.0/2.0)*b%b%k13+(1.0/2.0)*a%a%k31+a%b%k22;
hh2=k02+k03%b+(1.0/2.0)*k04%b%b +k12%a+a%b%k13+(1.0/2.0)*a%a%k22;
anew =a+(hh%gg2-gg%hh2)/(gg1%hh2-gg2%hh1);
bnew =b-(hh%gg1-gg%hh1)/(gg1%hh2-gg2%hh1);
abser =(pow(anew-a,2)+pow(bnew-b,2));
a=anew;
b=bnew;
ind=ind+1;
}
vec dens2=-log(2*3.141592653589793)-0.5*log(gg1%hh2-gg2%hh1)+(k10%a+k01%b+(1.0/2.0)*k20%a%a+(1.0/2.0)*k02%b%b+(1.0/6.0)*k30%a%a%a+(1.0/6.0)*k03%b%b%b+(1.0/24.0)*k40%a%a%a%a+(1.0/24.0)*k04%b%b%b%b+k11%a%b+(1.0/2.0)*k12%a%b%b+(1.0/2.0)*k21%a%a%b+(1.0/6.0)*a%b%b%b%k13+(1.0/6.0)*a%a%a%b%k31+(1.0/4.0)*a%a%b%b%k22-a%Xt-b%Yt);
abser=0.1+0*abser;
a.ones();
b.ones();
det=(kk10%kk01-kk11%kk11);
a=-(Xt-kk10)%kk20/det+(Yt-kk01)%kk11/det;
b=+(Xt-kk10)%kk11/det-(Yt-kk01)%kk02/det;
ind=0;
while((max(abser)>0.001)&&(ind<1500))
{
gg=kk10+kk20%a+(1.0/2.0)*kk30%a%a+(1.0/6.0)*kk40%a%a%a +kk11%b +(1.0/2.0)*kk12%b%b+kk21%a%b+(1.0/6.0)*b%b%b%kk13+(1.0/2.0)*a%a%b%kk31+(1.0/2.0)*a%b%b%kk22-Xt;
hh=kk01+kk02%b+(1.0/2.0)*kk03%b%b+(1.0/6.0)*kk04%b%b%b +kk11%a +kk12%a%b+(1.0/2.0)*kk21%a%a+(1.0/2.0)*a%b%b%kk13+(1.0/6.0)*a%a%a%kk31+(1.0/2.0)*a%a%b%kk22-Yt;
gg1=kk20+kk30%a+(1.0/2.0)*kk40%a%a+kk21%b+a%b%kk31+(1.0/2.0)*b%b%kk22;
gg2=kk11 +kk12%b+kk21%a+(1.0/2.0)*b%b%kk13+(1.0/2.0)*a%a%kk31+a%b%kk22;
hh1=kk11 +kk12%b+kk21%a+(1.0/2.0)*b%b%kk13+(1.0/2.0)*a%a%kk31+a%b%kk22;
hh2=kk02+kk03%b+(1.0/2.0)*kk04%b%b +kk12%a+a%b%kk13+(1.0/2.0)*a%a%kk22;
anew =a+(hh%gg2-gg%hh2)/(gg1%hh2-gg2%hh1);
bnew =b-(hh%gg1-gg%hh1)/(gg1%hh2-gg2%hh1);
abser =(pow(anew-a,2)+pow(bnew-b,2));
a=anew;
b=bnew;
ind=ind+1;
}
vec dens1=-log(2*3.141592653589793)-0.5*log(gg1%hh2-gg2%hh1)+(kk10%a+kk01%b+(1.0/2.0)*kk20%a%a+(1.0/2.0)*kk02%b%b+(1.0/6.0)*kk30%a%a%a+(1.0/6.0)*kk03%b%b%b+(1.0/24.0)*kk40%a%a%a%a+(1.0/24.0)*kk04%b%b%b%b+kk11%a%b+(1.0/2.0)*kk12%a%b%b+(1.0/2.0)*kk21%a%a%b+(1.0/6.0)*a%b%b%b%kk13+(1.0/6.0)*a%a%a%b%kk31+(1.0/4.0)*a%a%b%b%kk22-a%Xt-b%Yt);
resss.col(1)=a;
resss.col(2)=b;
resss.col(0)=log(exp(dens1)%probs+exp(dens2)%(1-probs));
List ret;
ret["like"] = resss;
ret["like2"] = dens1;
ret["like3"] = dens2;
ret["max"] = whch;
ret["probs"] = probs;
return(ret); }'
}
if(Dtype=='Edgeworth')
{
txtC='
vec x = (Xt-x0.col(0))/sqrt(x0.col(1));
vec y = (Yt-x0.col(4))/sqrt(x0.col(5));
vec rho30 =x0.col(2)/(pow(x0.col(1),1.5));
vec rho40 =x0.col(3)/(pow(x0.col(1),2));
vec rho03 =x0.col(6)/(pow(x0.col(5),1.5));
vec rho04 =x0.col(7)/(pow(x0.col(5),2));
vec rho11 =x0.col(8)/(sqrt(x0.col(1))%sqrt(x0.col(5)));
vec rho12 =x0.col(9)/(sqrt(x0.col(1))%x0.col(5));
vec rho21 =x0.col(10)/(x0.col(1)%sqrt(x0.col(5)));
vec rho22 =x0.col(11)/(x0.col(1)%x0.col(5));
vec rho13 =x0.col(12)/(sqrt(x0.col(1))%pow(x0.col(5),1.5));
vec rho31 =x0.col(13)/(pow(x0.col(1),1.5)%sqrt(x0.col(5)));
vec Hx1 = x;
vec Hx2 = pow(x,2)-1;
vec Hx3 = pow(x,3)-3*x;
vec Hx4 = pow(x,4)-6*pow(x,2)+3;
vec Hx5 = pow(x,5)-10*pow(x,3)+15*x;
vec Hx6 = pow(x,6)-15*pow(x,4)+45*pow(x,2)-15;
vec Hy1 = y;
vec Hy2 = pow(y,2)-1;
vec Hy3 = pow(y,3)-3*y;
vec Hy4 = pow(y,4)-6*pow(y,2)+3;
vec Hy5 = pow(y,5)-10*pow(y,3)+15*y;
vec Hy6 = pow(y,6)-15*pow(y,4)+45*pow(y,2)-15;
vec A=Hx3%rho30+3*Hx2%Hy1%rho21+3*Hx1%Hy2%rho12+Hy3%rho03;
vec B=Hx4%rho40+4*Hx3%Hy1%rho31+6*Hx2%Hy2%rho22+4*Hx1%Hy3%rho13+Hy4%rho04;
vec C=Hx6%pow(rho30,2)+6*Hx5%Hy1%rho21%rho30+6*Hx4%Hy2%rho12%rho30+2*Hx3%Hy3%rho03%rho30+9*Hx4%Hy2%pow(rho21,2)+18*Hx3%Hy3%rho12%rho21+6*Hx2%Hy4%rho03%rho21+9*Hx2%Hy4%pow(rho12,2)+6*Hx1%Hy5%rho03%rho12+Hy6%pow(rho03,2);
resss.col(0)=-log(2*pi*sqrt(x0.col(1)%x0.col(5))%sqrt(1 - pow(rho11,2)))+(-(pow((Xt-x0.col(0)),2)/x0.col(1) -2*rho11%(Xt-x0.col(0))%(Yt-x0.col(4))/sqrt(x0.col(1)%x0.col(5))+pow((Yt-x0.col(4)),2)/x0.col(5))/(2*(1-pow(rho11,2))))+log(abs(1+0.1666667*A+0.04166667*B+0.01388889*C));
List ret;
ret["like"] = resss;
ret["max"] = whch;
return(ret);
}'
}
if(Dtype=='Normal')
{
txtC='
vec probs=exp(-x0.col(28)-x0.col(29)-x0.col(30)-x0.col(31)-x0.col(32)-x0.col(33));
vec m00 = (1+0*x0.col(0));
vec m10 = (x0.col(0) -probs%x0.col(14))/(1-probs);
vec m20 = (x0.col(1) -probs%x0.col(15))/(1-probs);
vec m30 = (x0.col(2) -probs%x0.col(16))/(1-probs);
vec m40 = (x0.col(3) -probs%x0.col(17))/(1-probs);
vec m01 = (x0.col(4) -probs%x0.col(18))/(1-probs);
vec m02 = (x0.col(5) -probs%x0.col(19))/(1-probs);
vec m03 = (x0.col(6) -probs%x0.col(20))/(1-probs);
vec m04 = (x0.col(7) -probs%x0.col(21))/(1-probs);
vec m11 = (x0.col(8) -probs%x0.col(22))/(1-probs);
vec m12 = (x0.col(9) -probs%x0.col(23))/(1-probs);
vec m21 = (x0.col(10)-probs%x0.col(24))/(1-probs);
vec m22 = (x0.col(11)-probs%x0.col(25))/(1-probs);
vec m13 = (x0.col(12)-probs%x0.col(26))/(1-probs);
vec m31 = (x0.col(13)-probs%x0.col(27))/(1-probs);
vec mm00 = (1+0*x0.col(0));
vec mm10 = x0.col(14);
vec mm20 = x0.col(15);
vec mm30 = x0.col(16);
vec mm40 = x0.col(17);
vec mm01 = x0.col(18);
vec mm02 = x0.col(19);
vec mm03 = x0.col(20);
vec mm04 = x0.col(21);
vec mm11 = x0.col(22);
vec mm12 = x0.col(23);
vec mm21 = x0.col(24);
vec mm22 = x0.col(25);
vec mm13 = x0.col(26);
vec mm31 = x0.col(27);
vec k10 = m10;
vec k20 = m20-pow(m10,2);
vec k30 = m30-3*m20%m10+2*pow(m10,3);
vec k40 = m40 -4*m30%m10-3*pow(m20,2)+12*m20%pow(m10,2)-6*pow(m10,4);
vec k01 = m01;
vec k02 = m02- pow(m01,2) ;
vec k03 = m03-3*m02%m01+2*pow(m01,3) ;
vec k04 = m04-4*m03%m01-3*pow(m02,2)+12*m02%pow(m01,2)-6*pow(m01,4);
vec k11 = m11-m10%m01;
vec k21 = m21-2*m11%m10-m20%m01+2*pow(m10,2)%m01;
vec k12 = m12-2*m11%m01-m02%m10+2*pow(m01,2)%m10;
vec k22 = m22-2*m21%m01-2*m12%m10-m20%m02-2*pow(m11,2)+8*m11%m01%m10+2*m02%pow(m10,2)+2*m20%pow(m01,2)-6*pow(m10,2)%pow(m01,2) ;
vec k31 = m31-3*m21%m10-m30%m01-3*m20%m11+6*m11%pow(m10,2)+6*m20%m10%m01-6*pow(m10,3)%m01 ;
vec k13 = m13-3*m12%m01-m03%m10-3*m02%m11+6*m11%pow(m01,2)+6*m02%m01%m10-6*pow(m01,3)%m10;
vec kk10 = mm10;
vec kk20 = mm20-pow(mm10,2);
vec kk30 = mm30-3*mm20%mm10+2*pow(mm10,3);
vec kk40 = mm40 -4*mm30%mm10-3*pow(mm20,2)+12*mm20%pow(mm10,2)-6*pow(mm10,4);
vec kk01 = mm01;
vec kk02 = mm02- pow(mm01,2) ;
vec kk03 = mm03-3*mm02%mm01+2*pow(mm01,3) ;
vec kk04 = mm04-4*mm03%mm01-3*pow(mm02,2)+12*mm02%pow(mm01,2)-6*pow(mm01,4);
vec kk11 = mm11-mm10%mm01;
vec kk21 = mm21-2*mm11%mm10-mm20%mm01+2*pow(mm10,2)%mm01;
vec kk12 = mm12-2*mm11%mm01-mm02%mm10+2*pow(mm01,2)%mm10;
vec kk22 = mm22-2*mm21%mm01-2*mm12%mm10-mm20%mm02-2*pow(mm11,2)+8*mm11%mm01%mm10+2*mm02%pow(mm10,2)+2*mm20%pow(mm01,2)-6*pow(mm10,2)%pow(mm01,2) ;
vec kk31 = mm31-3*mm21%mm10-mm30%mm01-3*mm20%mm11+6*mm11%pow(mm10,2)+6*mm20%mm10%mm01-6*pow(mm10,3)%mm01 ;
vec kk13 = mm13-3*mm12%mm01-mm03%mm10-3*mm02%mm11+6*mm11%pow(mm01,2)+6*mm02%mm01%mm10-6*pow(mm01,3)%mm10;
vec det=(kk20%kk02-kk11%kk11);
vec dens1=-log(2*3.141592653589793)-0.5*log(abs(det))-0.5*((Xt-kk10)%(Xt-kk10)%kk02/det-(Xt-kk10)%(Yt-kk01)%kk11/det-(Xt-kk10)%(Yt-kk01)%kk11/det+(Yt-kk01)%(Yt-kk01)%kk20/det);
det=(k20%k02-k11%k11);
vec dens2=-log(2*3.141592653589793)-0.5*log(abs(det))-0.5*((Xt-k10)%(Xt-k10)%k02/det-(Xt-k10)%(Yt-k01)%k11/det-(Xt-k10)%(Yt-k01)%k11/det+(Yt-k01)%(Yt-k01)%k20/det);
resss.col(0)=log(exp(dens1)%probs+exp(dens2)%(1-probs));
List ret;
ret["like"] = resss;
ret["max"] = whch;
ret["probs"] = probs;
return(ret);
}'
}
}
#==============================================================================
# Generate TYPE of Solution
#==============================================================================
if(state4)
{
# DATA RESOLUTION -------------------------------------------------------------
if(!homo.res)
{
delt=cbind(diff(T.seq)/mesh,diff(T.seq)/mesh,diff(T.seq)/mesh,diff(T.seq)/mesh,diff(T.seq)/mesh,diff(T.seq)/mesh,diff(T.seq)/mesh,diff(T.seq)/mesh,
diff(T.seq)/mesh,diff(T.seq)/mesh,diff(T.seq)/mesh,diff(T.seq)/mesh,diff(T.seq)/mesh,diff(T.seq)/mesh,
diff(T.seq)/mesh,diff(T.seq)/mesh,diff(T.seq)/mesh,diff(T.seq)/mesh,diff(T.seq)/mesh,diff(T.seq)/mesh,diff(T.seq)/mesh,diff(T.seq)/mesh,
diff(T.seq)/mesh,diff(T.seq)/mesh,diff(T.seq)/mesh,diff(T.seq)/mesh,diff(T.seq)/mesh,diff(T.seq)/mesh,
diff(T.seq)/mesh,diff(T.seq)/mesh,diff(T.seq)/mesh,diff(T.seq)/mesh,diff(T.seq)/mesh,diff(T.seq)/mesh)
}
t=seq(0,100,1/100)
tro=34
if(Jdist=='Normal')
{
Mstar1 =
c("+1*mm1*m00"
,"+2*mm1*m10+1*mm2"
,"+3*mm1*m20+3*mm2*m10+1*mm3"
,"+4*mm1*m30+6*mm2*m20+4*mm3*m10+1*mm4"
,"+1*oo1*m00"
,"+2*oo1*m01+1*oo2"
,"+3*oo1*m02+3*oo2*m01+1*oo3"
,"+4*oo1*m03+6*oo2*m02+4*oo3*m01+1*oo4"
,"mm1*m01+oo1*m10+mm1*oo1"
,"mm1*m02+2*oo1*m11+2*mm1*oo1*m01+oo2*m10+mm1*oo2"
,"2*mm1*m11+mm2*m01+oo1*m20+2*mm1*oo1*m10+mm2*oo1"
,"2*mm1*m12+mm2*m02+2*oo1*m21+4*mm1*oo1*m11+2*mm2*oo1*m01+oo2*m20+2*mm1*oo2*m10+mm2*oo2"
,"mm1*m03+3*oo1*m12+3*mm1*oo1*m02+3*oo2*m11+3*mm1*oo2*m01+oo3*m10+mm1*oo3"
,"3*mm1*m21+3*mm2*m11+mm3*m01+oo1*m30+3*mm1*oo1*m20+3*mm2*oo1*m10+mm3*oo1")
Mstar2 =
c("+1*nn1*m00"
,"+2*nn1*m10+1*nn2"
,"+3*nn1*m20+3*nn2*m10+1*nn3"
,"+4*nn1*m30+6*nn2*m20+4*nn3*m10+1*nn4"
,"+1*pp1*m00"
,"+2*pp1*m01+1*pp2"
,"+3*pp1*m02+3*pp2*m01+1*pp3"
,"+4*pp1*m03+6*pp2*m02+4*pp3*m01+1*pp4"
,"nn1*m01+pp1*m10+nn1*pp1"
,"nn1*m02+2*pp1*m11+2*nn1*pp1*m01+pp2*m10+nn1*pp2"
,"2*nn1*m11+nn2*m01+pp1*m20+2*nn1*pp1*m10+nn2*pp1"
,"2*nn1*m12+nn2*m02+2*pp1*m21+4*nn1*pp1*m11+2*nn2*pp1*m01+pp2*m20+2*nn1*pp2*m10+nn2*pp2"
,"nn1*m03+3*pp1*m12+3*nn1*pp1*m02+3*pp2*m11+3*nn1*pp2*m01+pp3*m10+nn1*pp3"
,"3*nn1*m21+3*nn2*m11+nn3*m01+pp1*m30+3*nn1*pp1*m20+3*nn2*pp1*m10+nn3*pp1")
Mstar1x =
c("+1*mm1*m10"
,"+2*mm1*m20+1*mm2*m10"
,"+3*mm1*m30+3*mm2*m20+1*mm3*m10"
,"+4*mm1*m40+6*mm2*m30+4*mm3*m20+1*mm4*m10"
,"+1*oo1*m10"
,"+2*oo1*m11+1*oo2*m10"
,"+3*oo1*m12+3*oo2*m11+1*oo3*m10"
,"+4*oo1*m13+6*oo2*m12+4*oo3*m11+1*oo4*m10"
,"mm1*m11+oo1*m20+mm1*oo1*m10"
,"mm1*m12+2*oo1*m21+2*mm1*oo1*m11+oo2*m20+mm1*oo2*m10"
,"2*mm1*m21+mm2*m11+oo1*m30+2*mm1*oo1*m20+mm2*oo1*m10"
,"2*mm1*m22+mm2*m12+2*oo1*m31+4*mm1*oo1*m21+2*mm2*oo1*m11+oo2*m30+2*mm1*oo2*m20+mm2*oo2*m10"
,"mm1*m13+3*oo1*m22+3*mm1*oo1*m12+3*oo2*m21+3*mm1*oo2*m11+oo3*m20+mm1*oo3*m10"
,"3*mm1*m31+3*mm2*m21+mm3*m11+oo1*m40+3*mm1*oo1*m30+3*mm2*oo1*m20+mm3*oo1*m10")
Mstar1y =
c("+1*mm1*m01"
,"+2*mm1*m11+1*mm2*m01"
,"+3*mm1*m21+3*mm2*m11+1*mm3*m01"
,"+4*mm1*m31+6*mm2*m21+4*mm3*m11+1*mm4*m01"
,"+1*oo1*m01"
,"+2*oo1*m02+1*oo2*m01"
,"+3*oo1*m03+3*oo2*m02+1*oo3*m01"
,"+4*oo1*m04+6*oo2*m03+4*oo3*m02+1*oo4*m01"
,"mm1*m02+oo1*m11+mm1*oo1*m01"
,"mm1*m03+2*oo1*m11+2*mm1*oo1*m02+oo2*m11+mm1*oo2*m01"
,"2*mm1*m12+mm2*m02+oo1*m21+2*mm1*oo1*m11+mm2*oo1*m01"
,"2*mm1*m13+mm2*m03+2*oo1*m22+4*mm1*oo1*m12+2*mm2*oo1*m02+oo2*m21+2*mm1*oo2*m11+mm2*oo2*m01"
,"mm1*m04+3*oo1*m13+3*mm1*oo1*m03+3*oo2*m11+3*mm1*oo2*m02+oo3*m11+mm1*oo3*m01"
,"3*mm1*m22+3*mm2*m12+mm3*m02+oo1*m31+3*mm1*oo1*m21+3*mm2*oo1*m11+mm3*oo1*m02")
Mstar2x =
c("+1*nn1*m10"
,"+2*nn1*m20+1*nn2*m10"
,"+3*nn1*m30+3*nn2*m20+1*nn3*m10"
,"+4*nn1*m40+6*nn2*m30+4*nn3*m20+1*nn4*m10"
,"+1*pp1*m10"
,"+2*pp1*m11+1*pp2*m10"
,"+3*pp1*m12+3*pp2*m11+1*pp3*m10"
,"+4*pp1*m13+6*pp2*m12+4*pp3*m11+1*pp4*m10"
,"nn1*m11+pp1*m20+nn1*pp1*m10"
,"nn1*m12+2*pp1*m21+2*nn1*pp1*m11+pp2*m20+nn1*pp2*m10"
,"2*nn1*m21+nn2*m11+pp1*m30+2*nn1*pp1*m20+nn2*pp1*m10"
,"2*nn1*m22+nn2*m12+2*pp1*m31+4*nn1*pp1*m21+2*nn2*pp1*m11+pp2*m30+2*nn1*pp2*m20+nn2*pp2*m10"
,"nn1*m13+3*pp1*m22+3*nn1*pp1*m12+3*pp2*m21+3*nn1*pp2*m11+pp3*m20+nn1*pp3*m10"
,"3*nn1*m31+3*nn2*m21+nn3*m11+pp1*m40+3*nn1*pp1*m30+3*nn2*pp1*m20+nn3*pp1*m10")
Mstar2y =
c("+1*nn1*m01"
,"+2*nn1*m11+1*nn2*m01"
,"+3*nn1*m21+3*nn2*m11+1*nn3*m01"
,"+4*nn1*m31+6*nn2*m21+4*nn3*m11+1*nn4*m01"
,"+1*pp1*m01"
,"+2*pp1*m02+1*pp2*m01"
,"+3*pp1*m03+3*pp2*m02+1*pp3*m01"
,"+4*pp1*m04+6*pp2*m03+4*pp3*m02+1*pp4*m01"
,"nn1*m02+pp1*m11+nn1*pp1*m01"
,"nn1*m03+2*pp1*m12+2*nn1*pp1*m02+pp2*m11+nn1*pp2*m01"
,"2*nn1*m12+nn2*m02+pp1*m21+2*nn1*pp1*m11+nn2*pp1*m01"
,"2*nn1*m13+nn2*m03+2*pp1*m22+4*nn1*pp1*m12+2*nn2*pp1*m02+pp2*m21+2*nn1*pp2*m11+nn2*pp2*m01"
,"nn1*m04+3*pp1*m13+3*nn1*pp1*m03+3*pp2*m12+3*nn1*pp2*m02+pp3*m11+nn1*pp3*m01"
,"3*nn1*m22+3*nn2*m12+nn3*m02+pp1*m31+3*nn1*pp1*m21+3*nn2*pp1*m11+nn3*pp1*m01")
}
if(Jdist=='MVNormal')
{
if(Jtype=='Add')
{
Mstar1 =
c("+1*mv10*m00"
,"+2*mv10*m10+1*mv20"
,"+3*mv10*m20+3*mv20*m10+1*mv30"
,"+4*mv10*m30+6*mv20*m20+4*mv30*m10+1*mv40"
,"+1*mv01*m00"
,"+2*mv01*m01+1*mv02"
,"+3*mv01*m02+3*mv02*m01+1*mv03"
,"+4*mv01*m03+6*mv02*m02+4*mv03*m01+1*mv04"
,"mv10*m01+mv01*m10+mv11"
,"mv10*m02+2*mv01*m11+2*mv11*m01+mv02*m10+mv12"
,"2*mv10*m11+mv20*m01+mv01*m20+2*mv11*m10+mv21"
,"2*mv10*m12+mv20*m02+2*mv01*m21+4*mv11*m11+2*mv21*m01+mv02*m20+2*mv12*m10+mv22"
,"mv10*m03+3*mv01*m12+3*mv11*m02+3*mv02*m11+3*mv12*m01+mv03*m10+mv13"
,"3*mv10*m21+3*mv20*m11+mv30*m01+mv01*m30+3*mv11*m20+3*mv21*m10+mv31")
}
if(Jtype=='Mult')
{
Mstar1 =
c("m10*(mv10)"
,"m20*(mv20+2*mv10)"
,"m30*(mv30+3*mv20+3*mv10)"
,"m40*(mv40+4*mv30+6*mv20+4*mv10)"
,"m01*(mv01)"
,"m02*(mv02+2*mv01)"
,"m03*(mv03+3*mv02+3*mv01)"
,"m04*(mv04+4*mv03+6*mv02+4*mv01)"
,"m11*(mv11+mv01+mv10)"
,"m12*(mv12+mv02+2*mv11+2*mv01+mv10)"
,"m21*(mv21+2*mv11+mv01+mv20+2*mv10)"
,"m22*(mv22+2*mv12+mv02+2*mv21+4*mv11+2*mv01+mv20+2*mv10)"
,"m13*(mv13+mv03+3*mv12+3*mv02+3*mv11+3*mv01+mv10)"
,"m31*(mv31+3*mv21+3*mv11+mv01+mv30+3*mv20+3*mv10)")
}
}
MAT=rbind(
c("1*m00","1*m10","1*m20","1*m01","1*m02","1*m11","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","")
,c("2*m10","2*m20","2*m30","2*m11","2*m12","2*m21","","","","","","","1*m00","1*m10","1*m20","1*m01","1*m02","1*m11","","","","","","","","","","","","","","","","","","")
,c("3*m20","3*m30","3*m40","3*m21","3*m22","3*m31","","","","","","","3*m10","3*m20","3*m30","3*m11","3*m12","3*m21","","","","","","","","","","","","","","","","","","")
,c("4*m30","4*m40","4*m50","4*m31","4*m32","4*m41","","","","","","","6*m20","6*m30","6*m40","6*m21","6*m22","6*m31","","","","","","","","","","","","","","","","","","")
,c("","","","","","","1*m00","1*m10","1*m20","1*m01","1*m02","1*m11","","","","","","","","","","","","","","","","","","","","","","","","")
,c("","","","","","","2*m01","2*m11","2*m21","2*m02","2*m03","2*m12","","","","","","","","","","","","","","","","","","","1*m00","1*m10","1*m20","1*m01","1*m02","1*m11")
,c("","","","","","","3*m02","3*m12","3*m22","3*m03","3*m04","3*m13","","","","","","","","","","","","","","","","","","","3*m01","3*m11","3*m21","3*m02","3*m03","3*m12")
,c("","","","","","","4*m03","4*m13","4*m23","4*m04","4*m05","4*m14","","","","","","","","","","","","","","","","","","","6*m02","6*m12","6*m22","6*m03","6*m04","6*m13")
,c("1*m01","1*m11","1*m21","1*m02","1*m03","1*m12","1*m10","1*m20","1*m30","1*m11","1*m12","1*m21","","","","","","","1*m00","1*m10","1*m20","1*m01","1*m02","1*m11","1*m00","1*m10","1*m20","1*m01","1*m02","1*m11","","","","","","")
,c("1*m02","1*m12","1*m22","1*m03","1*m04","1*m13","2*m11","2*m21","2*m31","2*m12","2*m13","2*m22","","","","","","","2*m01","2*m11","2*m21","2*m02","2*m03","2*m12","2*m01","2*m11","2*m21","2*m02","2*m03","2*m12","1*m10","1*m20","1*m30","1*m11","1*m12","1*m21")
,c("2*m11","2*m21","2*m31","2*m12","2*m13","2*m22","1*m20","1*m30","1*m40","1*m21","1*m22","1*m31","1*m01","1*m11","1*m21","1*m02","1*m03","1*m12","2*m10","2*m20","2*m30","2*m11","2*m12","2*m21","2*m10","2*m20","2*m30","2*m11","2*m12","2*m21","","","","","","")
,c("2*m12","2*m22","2*m32","2*m13","2*m14","2*m23","2*m21","2*m31","2*m41","2*m22","2*m23","2*m32","1*m02","1*m12","1*m22","1*m03","1*m04","1*m13","4*m11","4*m21","4*m31","4*m12","4*m13","4*m22","4*m11","4*m21","4*m31","4*m12","4*m13","4*m22","1*m20","1*m30","1*m40","1*m21","1*m22","1*m31")
,c("1*m03","1*m13","1*m23","1*m04","1*m05","1*m14","3*m12","3*m22","3*m32","3*m13","3*m14","3*m23","","","","","","","3*m02","3*m12","3*m22","3*m03","3*m04","3*m13","3*m02","3*m12","3*m22","3*m03","3*m04","3*m13","3*m11","3*m21","3*m31","3*m12","3*m13","3*m22")
,c("3*m21","3*m31","3*m41","3*m22","3*m23","3*m32","1*m30","1*m40","1*m50","1*m31","1*m32","1*m41","3*m11","3*m21","3*m31","3*m12","3*m13","3*m22","3*m20","3*m30","3*m40","3*m21","3*m22","3*m31","3*m20","3*m30","3*m40","3*m21","3*m22","3*m31","","","","","","")
)
MAT2=rbind(
c("1*mm00","1*mm10","1*mm20","1*mm01","1*mm02","1*mm11","","","","","","","","","","","","","","","","","","","","","","","","","","","","","","")
,c("2*mm10","2*mm20","2*mm30","2*mm11","2*mm12","2*mm21","","","","","","","1*mm00","1*mm10","1*mm20","1*mm01","1*mm02","1*mm11","","","","","","","","","","","","","","","","","","")
,c("3*mm20","3*mm30","3*mm40","3*mm21","3*mm22","3*mm31","","","","","","","3*mm10","3*mm20","3*mm30","3*mm11","3*mm12","3*mm21","","","","","","","","","","","","","","","","","","")
,c("4*mm30","4*mm40","4*mm50","4*mm31","4*mm32","4*mm41","","","","","","","6*mm20","6*mm30","6*mm40","6*mm21","6*mm22","6*mm31","","","","","","","","","","","","","","","","","","")
,c("","","","","","","1*mm00","1*mm10","1*mm20","1*mm01","1*mm02","1*mm11","","","","","","","","","","","","","","","","","","","","","","","","")
,c("","","","","","","2*mm01","2*mm11","2*mm21","2*mm02","2*mm03","2*mm12","","","","","","","","","","","","","","","","","","","1*mm00","1*mm10","1*mm20","1*mm01","1*mm02","1*mm11")
,c("","","","","","","3*mm02","3*mm12","3*mm22","3*mm03","3*mm04","3*mm13","","","","","","","","","","","","","","","","","","","3*mm01","3*mm11","3*mm21","3*mm02","3*mm03","3*mm12")
,c("","","","","","","4*mm03","4*mm13","4*mm23","4*mm04","4*mm05","4*mm14","","","","","","","","","","","","","","","","","","","6*mm02","6*mm12","6*mm22","6*mm03","6*mm04","6*mm13")
,c("1*mm01","1*mm11","1*mm21","1*mm02","1*mm03","1*mm12","1*mm10","1*mm20","1*mm30","1*mm11","1*mm12","1*mm21","","","","","","","1*mm00","1*mm10","1*mm20","1*mm01","1*mm02","1*mm11","1*mm00","1*mm10","1*mm20","1*mm01","1*mm02","1*mm11","","","","","","")
,c("1*mm02","1*mm12","1*mm22","1*mm03","1*mm04","1*mm13","2*mm11","2*mm21","2*mm31","2*mm12","2*mm13","2*mm22","","","","","","","2*mm01","2*mm11","2*mm21","2*mm02","2*mm03","2*mm12","2*mm01","2*mm11","2*mm21","2*mm02","2*mm03","2*mm12","1*mm10","1*mm20","1*mm30","1*mm11","1*mm12","1*mm21")
,c("2*mm11","2*mm21","2*mm31","2*mm12","2*mm13","2*mm22","1*mm20","1*mm30","1*mm40","1*mm21","1*mm22","1*mm31","1*mm01","1*mm11","1*mm21","1*mm02","1*mm03","1*mm12","2*mm10","2*mm20","2*mm30","2*mm11","2*mm12","2*mm21","2*mm10","2*mm20","2*mm30","2*mm11","2*mm12","2*mm21","","","","","","")
,c("2*mm12","2*mm22","2*mm32","2*mm13","2*mm14","2*mm23","2*mm21","2*mm31","2*mm41","2*mm22","2*mm23","2*mm32","1*mm02","1*mm12","1*mm22","1*mm03","1*mm04","1*mm13","4*mm11","4*mm21","4*mm31","4*mm12","4*mm13","4*mm22","4*mm11","4*mm21","4*mm31","4*mm12","4*mm13","4*mm22","1*mm20","1*mm30","1*mm40","1*mm21","1*mm22","1*mm31")
,c("1*mm03","1*mm13","1*mm23","1*mm04","1*mm05","1*mm14","3*mm12","3*mm22","3*mm32","3*mm13","3*mm14","3*mm23","","","","","","","3*mm02","3*mm12","3*mm22","3*mm03","3*mm04","3*mm13","3*mm02","3*mm12","3*mm22","3*mm03","3*mm04","3*mm13","3*mm11","3*mm21","3*mm31","3*mm12","3*mm13","3*mm22")
,c("3*mm21","3*mm31","3*mm41","3*mm22","3*mm23","3*mm32","1*mm30","1*mm40","1*mm50","1*mm31","1*mm32","1*mm41","3*mm11","3*mm21","3*mm31","3*mm12","3*mm13","3*mm22","3*mm20","3*mm30","3*mm40","3*mm21","3*mm22","3*mm31","3*mm20","3*mm30","3*mm40","3*mm21","3*mm22","3*mm31","","","","","","")
)
namess = c('a00','a10','a20','a01','a02','a11',
'b00','b10','b20','b01','b02','b11',
'c00','c10','c20','c01','c02','c11',
'd00','d10','d20','d01','d02','d11',
'e00','e10','e20','e01','e02','e11',
'f00','f10','f20','f01','f02','f11')
objlist=objects(pos=1)
checknames = rep(0,length(namess))
for(i in 1:length(namess))
{
if(sum(objlist==namess[i])==1){checknames[i] = 1}
}
whichnames = which(checknames==1)
func.list.timehomo=rep(0,length(namess))
for(i in whichnames)
{
# which expressions vary over time
result=eval(body(namess[i]))
func.list.timehomo[i]=2-(sum(diff(result)==0)==(length(result)-1))
}
if(any(func.list.timehomo==2)){homo=F}
dims= rep('',14*2)
for(i in 1:14)
{
for(j in whichnames)
{
if(MAT[i,j]!='')
{
dims[i] =paste0(dims[i],'+(',body(namess[j])[2],')',c('*','%')[func.list.timehomo[j]],'(',MAT[i,j],')')
dims[i+14] =paste0(dims[i+14],'+(',body(namess[j])[2],')',c('*','%')[func.list.timehomo[j]],'(',MAT2[i,j],')')
}
}
}
#namess2 = c('Lam00','Lamy0','Nmu11','Nsig11','Nmu21','Nsig21','Nmu12','Nsig12','Nmu22','Nsig22','Lam10','Lamx2','Lamy1','Lamy2','Jmu1','Jmu2','Jsig11','Jsig12','Jsig22')
namess2 = c('Lam00','Lam10','Lam01','Jmu1','Jmu2','Jsig11','Jsig12','Jsig22')
checknames2 = rep(0,length(namess2))
for(i in 1:length(namess2))
{
if(sum(objlist==namess2[i])==1){checknames2[i] = 1}
}
whichnames2 = which(checknames2==1)
func.list.timehomo2=rep(0,length(namess2))
for(i in whichnames2)
{
# which expressions vary over time
result=eval(body(namess2[i]))
func.list.timehomo2[i]=2-(sum(diff(result)==0)==(length(result)-1))
}
if(any(func.list.timehomo2==2)){homo=F}
#print(func.list.timehomo)
#print(func.list.timehomo2)
if(checknames2[1]==0){Lam00=function(t){0}}
if(checknames2[2]==0){Lam10=function(t){0}}
if(checknames2[3]==0){Lam01=function(t){0}}
#if(checknames2[4]==0){Nsig11=function(t){0};Nmu11=function(t){0}}
#if(checknames2[5]==0){Nmu21=function(t){0}}
#if(checknames2[6]==0){Nsig21=function(t){0};Nmu21=function(t){0}}
#if(checknames2[7]==0){Nmu12=function(t){0}}
#if(checknames2[8]==0){Nsig12=function(t){0};Nmu12=function(t){0}}
#if(checknames2[9]==0){Nmu22=function(t){0}}
#if(checknames2[10]==0){Nsig22=function(t){0};Nmu22=function(t){0}}
#if(checknames2[11]==0){Lam10=function(t){0}}
#if(checknames2[12]==0){Lamx2=function(t){0}}
#if(checknames2[13]==0){Lamy1=function(t){0}}
#if(checknames2[14]==0){Lamy2=function(t){0}}
if(checknames2[4]==0){Jmu1=function(t){0}}
if(checknames2[5]==0){Jmu2=function(t){0}}
if(checknames2[6]==0){Jsig11=function(t){0}}
if(checknames2[7]==0){Jsig12=function(t){0}}
if(checknames2[8]==0){Jsig22=function(t){0}}
if(checknames2[1]==1)
{
for(i in 1:14)
{
dims[i] = paste0(dims[i],'+(',body(namess2[1])[2],')',c('*','%')[func.list.timehomo2[1]],'(',Mstar1[i],')')
}
}
if(checknames2[2]==1)
{
for(i in 1:14)
{
dims[i] = paste0(dims[i],'+(',body(namess2[2])[2],')',c('*','%')[func.list.timehomo2[2]],'(',Mstar2[i],')')
}
}
if(checknames2[3]==1)
{
for(i in 1:14)
{
dims[i] = paste0(dims[i],'+(',body(namess2[3])[2],')',c('*','%')[func.list.timehomo2[3]],'(',Mstar3[i],')')
}
}
# if(checknames2[11]==1)
#{
# for(i in 1:14)
# {
# dims[i] = paste0(dims[i],'+(',body(namess2[11])[2],')',c('*','%')[func.list.timehomo2[11]],'(',Mstar1x[i],')')
# }
#}
# if(checknames2[12]==1)
#{
# for(i in 1:14)
# {
# dims[i] = paste0(dims[i],'+(',body(namess2[12])[2],')',c('*','%')[func.list.timehomo2[12]],'(',Mstar2x[i],')')
# }
#}
# if(checknames2[13]==1)
#{
# for(i in 1:14)
# {
# dims[i] = paste0(dims[i],'+(',body(namess2[13])[2],')',c('*','%')[func.list.timehomo2[13]],'(',Mstar1y[i],')')
# }
#}
# if(checknames2[14]==1)
#{
# for(i in 1:14)
# {
# dims[i] = paste0(dims[i],'+(',body(namess2[14])[2],')',c('*','%')[func.list.timehomo2[14]],'(',Mstar2y[i],')')
# }
#}
#if((Jdist=='Normal'))
#{
#prem =
#'
# double mm1 = mu11 ;
# double mm2 = pow(mu11,2)+ pow(sig11,2) ;
# double mm3 = pow(mu11,3)+ 3* pow(mu11,1)*pow(sig11,2) ;
# double mm4 = pow(mu11,4)+ 6* pow(mu11,2)*pow(sig11,2)+3*pow(sig11,4) ;
#
# double oo1 = mu21 ;
# double oo2 = pow(mu21,2)+ pow(sig21,2) ;
# double oo3 = pow(mu21,3)+ 3* pow(mu21,1)*pow(sig21,2) ;
# double oo4 = pow(mu21,4)+ 6* pow(mu21,2)*pow(sig21,2)+3*pow(sig21,4) ;
#
# double nn1 = mu12 ;
# double nn2 = pow(mu12,2)+ pow(sig12,2) ;
# double nn3 = pow(mu12,3)+ 3* pow(mu12,1)*pow(sig12,2) ;
# double nn4 = pow(mu12,4)+ 6* pow(mu12,2)*pow(sig12,2)+3*pow(sig12,4) ;
# double pp1 = mu22 ;
# double pp2 = pow(mu22,2)+ pow(sig22,2) ;
# double pp3 = pow(mu22,3)+ 3* pow(mu22,1)*pow(sig22,2) ;
# double pp4 = pow(mu22,4)+ 6* pow(mu22,2)*pow(sig22,2)+3*pow(sig22,4) ;
#'
#}
if(Jdist=='MVNormal')
{
#print('MVNormal')
prem =
'
double mv10 = mu1 ;
double mv20 = pow(mu1,2)+ sig11 ;
double mv30 = pow(mu1,3)+ 3* pow(mu1,1)*sig11 ;
double mv40 = pow(mu1,4)+ 6* pow(mu1,2)*sig11+3*pow(sig11,2) ;
double mv01 = mu2 ;
double mv02 = pow(mu2,2)+ sig22 ;
double mv03 = pow(mu2,3)+ 3* pow(mu2,1)*sig22 ;
double mv04 = pow(mu2,4)+ 6* pow(mu2,2)*sig22+3*pow(sig22,2) ;
double mv11 = mu1*mu2+sig12 ;
double mv12 = mu1*pow(mu2,2)+2*mu2*sig12+mu1*sig22 ;
double mv21 = pow(mu1,2)*mu2+2*mu1*sig12+mu2*sig11 ;
double mv22 = pow(mu1,2)*pow(mu2,2)+pow(mu2,2)*sig11+pow(mu1,2)*sig22+4*mu1*mu2*sig12 +sig11*sig22+2*sig12*sig12;
double mv13 = mu1*pow(mu2,3)+3*pow(mu2,2)*sig12+3*mu1*mu2*sig22 +3*sig12*sig22;
double mv31 = mu2*pow(mu1,3)+3*pow(mu1,2)*sig12+3*mu1*mu2*sig11 +3*sig12*sig11;
'
}
premprem =rep('',4)
#'Lam00','Lamy0','Nmu11','Nsig11','Nmu21','Nsig21','Nmu12','Nsig12','Nmu22','Nsig22','Lam10','Lamx2','Lamy1','Lamy2'
#if(checknames2[3]==1){premprem[1]= paste0('double mu11=',body('Nmu11')[2],';','double sig11=',body('Nsig11')[2],';')}
#if(checknames2[5]==1){premprem[3]= paste0('double mu21=',body('Nmu21')[2],';','double sig21=',body('Nsig21')[2],';')}
#if(checknames2[7]==1){premprem[2]= paste0('double mu12=',body('Nmu12')[2],';','double sig12=',body('Nsig12')[2],';')}
#if(checknames2[9]==1){premprem[4]= paste0('double mu22=',body('Nmu22')[2],';','double sig22=',body('Nsig22')[2],';')}
#namess2 = c('Lam00','Lam10','Lam01','Jmu1','Jmu2','Jsig11','Jsig12','Jsig22')
#if(checknames2[3]==0){premprem[1]= paste0('double mu11=',0,';','double sig11=',0,';')}
#if(checknames2[5]==0){premprem[3]= paste0('double mu21=',0,';','double sig21=',0,';')}
#if(checknames2[7]==0){premprem[2]= paste0('double mu12=',0,';','double sig12=',0,';')}
#if(checknames2[9]==0){premprem[4]= paste0('double mu22=',0,';','double sig22=',0,';')}
if(Jdist=='MVNormal')
{
#print('MVNormal')
#namess2 = c('Lam00','Lam10','Lam01','Jmu1','Jmu2','Jsig11','Jsig12','Jsig22')
premprem[1]= paste0('double mu1=',body('Jmu1')[2],';\n','double mu2=',body('Jmu2')[2],';\n','double sig11=',body('Jsig11')[2],';\n','double sig12=',body('Jsig12')[2],';\n','double sig22=',body('Jsig22')[2],';\n')
premprem[2:4]=''
}
prm=''
for(i in 1:4)
{
prm = paste0(prm,premprem[i],'\n ')
}
b = rep('(+0*a.col(0))',6)
if(checknames2[1]==1){b[1] = paste0('(',body('Lam00')[2],')',c('*','%')[func.list.timehomo2[1 ]],'m00')}else{b[1] = paste0('(+0*a.col(0))')}
if(checknames2[2]==1){b[2] = paste0('(',body('Lam10')[2],')',c('*','%')[func.list.timehomo2[2 ]],'m10')}else{b[2] = paste0('(+0*a.col(0))')}
if(checknames2[3]==1){b[3] = paste0('(',body('Lam01')[2],')',c('*','%')[func.list.timehomo2[3 ]],'m01')}else{b[3] = paste0('(+0*a.col(0))')}
#if(checknames2[12]==1){b[4] = paste0('(',body('Lamx2')[2],')',c('*','%')[func.list.timehomo2[12]],'mm10')}else{b[4] = paste0('(+0*a.col(0))')}
#if(checknames2[13]==1){b[5] = paste0('(',body('Lamy1')[2],')',c('*','%')[func.list.timehomo2[13]],'mm01')}else{b[5] = paste0('(+0*a.col(0))')}
#if(checknames2[14]==1){b[6] = paste0('(',body('Lamy2')[2],')',c('*','%')[func.list.timehomo2[14]],'mm01')}else{b[6] = paste0('(+0*a.col(0))')}
for(i in 1:28)
{
dims[i]= paste0( 'atemp.col(',i-1,')=',dims[i],';')
}
for(i in 1:6)
{
b[i]= paste0( 'atemp.col(',28+i-1,')=',b[i],';')
}
pr="
vec m00 = (1+0*a.col(0));
vec m10 = a.col(0);
vec m20 = a.col(1);
vec m30 = a.col(2);
vec m40 = a.col(3);
vec m01 = a.col(4);
vec m02 = a.col(5);
vec m03 = a.col(6);
vec m04 = a.col(7);
vec m11 = a.col(8);
vec m12 = a.col(9);
vec m21 = a.col(10);
vec m22 = a.col(11);
vec m13 = a.col(12);
vec m31 = a.col(13);
vec mm00 = (1+0*a.col(0));
vec mm10 = a.col(14);
vec mm20 = a.col(15);
vec mm30 = a.col(16);
vec mm40 = a.col(17);
vec mm01 = a.col(18);
vec mm02 = a.col(19);
vec mm03 = a.col(20);
vec mm04 = a.col(21);
vec mm11 = a.col(22);
vec mm12 = a.col(23);
vec mm21 = a.col(24);
vec mm22 = a.col(25);
vec mm13 = a.col(26);
vec mm31 = a.col(27);
"
pr2="
vec m00 = (1+0*x0.col(0));
vec m10 = x0.col(0);
vec m20 = x0.col(1);
vec m30 = x0.col(2);
vec m40 = x0.col(3);
vec m01 = x0.col(4);
vec m02 = x0.col(5);
vec m03 = x0.col(6);
vec m04 = x0.col(7);
vec m11 = x0.col(8);
vec m12 = x0.col(9);
vec m21 = x0.col(10);
vec m22 = x0.col(11);
vec m13 = x0.col(12);
vec m31 = x0.col(13);
vec mm00 = (1+0*x0.col(0));
vec mm01 = x0.col(14);
vec mm02 = x0.col(15);
vec mm03 = x0.col(16);
vec mm04 = x0.col(17);
vec mm10 = x0.col(18);
vec mm20 = x0.col(19);
vec mm30 = x0.col(20);
vec mm40 = x0.col(21);
vec mm11 = x0.col(22);
vec mm12 = x0.col(23);
vec mm21 = x0.col(24);
vec mm22 = x0.col(25);
vec mm13 = x0.col(26);
vec mm31 = x0.col(27);
"
odekernel=paste0(pr,prm,prem,'\n',dims[1],'\n',dims[2],'\n',dims[3],'\n',dims[4],
'\n',dims[5],'\n',dims[6],'\n',dims[7],'\n',dims[8],'\n',dims[9],'\n',dims[10],
'\n',dims[11],'\n',dims[12],'\n',dims[13],'\n',dims[14],'\n',dims[1+14],'\n',dims[2+14],'\n',dims[3+14],'\n',dims[4+14],
'\n',dims[5+14],'\n',dims[6+14],'\n',dims[7+14],'\n',dims[8+14],'\n',dims[9+14],'\n',dims[10+14],
'\n',dims[11+14],'\n',dims[12+14],'\n',dims[13+14],'\n',dims[14+14],'\n',b[1],'\n',b[2],'\n',b[3],'\n',b[4],'\n',b[5],'\n',b[6])
# write.table(odekernel,'Check_ODE_dims.txt')
# WRIGHT AND SOURCE -----------------------------------------------------------
txt.full=paste(txtA,odekernel,txtB,txtC)
type.sol =" GENERALIZED QUADRATIC DIFFUSON"
}
#library(Rcpp)
#library(RcppArmadillo)
if(wrt)
{
write(txt.full,'BiJGQD.mcmc.cpp')
}
stre="Compiling C++ code. Please wait."
cat(stre, " \r")
flush.console()
sourceCpp(code=txt.full)
cat(' ','\r')
seq1=c(1,2,3,4,0,0,0,0,1,1,2,2,1,3)
seq2=c(0,0,0,0,1,2,3,4,1,2,1,2,3,1)
#==============================================================================
# Interface Module
#==============================================================================
trim <- function (x) gsub("([[:space:]])", "", x)
namess4=c('a00','a10','a20','a01','a02','a11',
'b00','b10','b20','b01','b02','b11',
'c00','c10','c20','c01','c02','c11',
'd00','d10','d20','d01','d02','d11',
'e00','e10','e20','e01','e02','e11',
'f00','f10','f20','f01','f02','f11')
indnames =rep(0,36)
for(i in 1:36)
{
if(sum(obs==namess4[i]))
{
indnames[i]=TRUE
namess4[i]=paste0(namess4[i],' : ',trim(body(namess4[i])[2]))
}
}
indnames2 = rep(0,5)
for(j in whichnames2)
{
if(sum(obs==namess2[j]))
{
indnames2[j]=TRUE
namess2[j]=paste0(namess2[j],' : ',trim(body(namess2[j])[2]))
}
}
prior.list = trim(prior.list)
namess4=matrix(namess4,length(namess4),1)
buffer0=c('================================================================')
buffer1=c('----------------------------------------------------------------')
buffer2=c('................................................................')
buffer3=c('... ... ... ... ... ... ... ... ... ... ... ')
buffer4=c('_____________________ Drift Coefficients _______________________')
buffer5=c('___________________ Diffusion Coefficients _____________________')
buffer6=c('_____________________ Prior Distributions ______________________')
buffer7=c('_______________________ Model/Chain Info _______________________')
buffer8=c('......................... Intensity ............................')
buffer9=c('........................... Jumps ..............................')
buffer10=c('_______________________ Jump Components ________________________')
Info=c(buffer0,type.sol,buffer0,buffer4,
namess4[1:6][which(indnames[1:6]==T)],
buffer3,
namess4[7:12][which(indnames[7:12]==T)],
buffer5,
namess4[13:18][which(indnames[13:18]==T)],
buffer3,
namess4[19:24][which(indnames[19:24]==T)],
buffer3,
namess4[25:30][which(indnames[25:30]==T)],
buffer3,
namess4[31:36][which(indnames[31:36]==T)],
buffer10,
buffer8,
namess2[c(1:3)][which(indnames2[1:3]==T)],
buffer9,
namess2[c(4:8)][which(indnames2[4:8]==T)],
buffer6,'',prior.list)
Info=data.frame(matrix(Info,length(Info),1))
colnames(Info)=''
if(print.output)
{
print(Info,row.names = FALSE,right=F)
}
############################################################################
############################################################################
############################################################################
tme=Sys.time()
par.matrix=matrix(0,length(theta),updates)
probs=rep(0,updates)
ll=rep(0,updates)
errors =ll
acc=ll
kk=0
par.matrix[,1]=theta
prop.matrix =par.matrix
retries = 0
max.retries = 0
retry.count = 0
retry.indexes = c()
success = TRUE
strts =matrix(0,nnn-1,2)
strts[,1]=rtf[1]
strts[,2]=rtf[2]
rs=solver(X1[-nnn],X2[-nnn],X1[-1],X2[-1],c(0,theta),mesh,delt,nnn-1,T.seq[-nnn],strts,tro,secmom,seq1,seq2)
if(is.na(sum(rs$like[,1])))
{
retry.count =1
while(is.na(sum(rs$like[,1]))&&(retry.count<=10))
{
theta.start=theta+rnorm(length(theta),sd=sds)
rs = solver(X1[-nnn],X2[-nnn],X1[-1],X2[-1],c(0,theta.start),mesh,delt,nnn-1,T.seq[-nnn],strts,tro,secmom,seq1,seq2)
if(is.na(sum(rs$like[,1]))){retry.count=retry.count+1}
}
}
lold=sum(rs$like[,1])
# print(lold)
errors[1] = rs$max
if(is.na(lold)){print('Fail: Could not evaluate likelihood at initial perameters.');failed.chain=T;}
ll[1]=lold
excess = matrix(0,2,updates)
term.discont = rep(0,nnn-1)
decodes = term.discont
decodes.prob = decodes
dec.matrix =matrix(0,nnn-1,updates)
ltemp1 = rs$like3+log(1-rs$probs)
ltemp2 = rs$like2+log(rs$probs)
muvec = theta
covvec = diag(1/(adapt^2/length(theta))*(sds)^2)
if(adapt==0)
{
pb <- txtProgressBar(1,updates,1,style = 1,width = 65)
failed.chain=F
i=2
while(i<=updates)
{
theta.temp=theta
theta=theta+rnorm(length(theta),sd=sds)
prop.matrix[,i] = theta
tempp=solver(X1[-nnn],X2[-nnn],X1[-1],X2[-1],c(0,theta),mesh,delt,nnn-1,T.seq[-nnn],strts,tro,secmom,seq1,seq2)
strts=tempp$like[,2:3]*recycle+strts*(1-recycle)
errors[i] = tempp$max
lnew=sum(tempp$like[,1])
rat=min(exp(lnew-lold)*pp(theta)/pp(theta.temp),1)
if(is.na(rat))
{
retry.count =1
retries=retries+1
retry.indexes[retries] = i
max.retries=max.retries+1
while(is.na(rat)&&(retry.count<=10))
{
i = max(i-10,2)
theta = par.matrix[,i]
theta=theta+rnorm(length(theta),sd=sds)
prop.matrix[,i] = theta
tempp=solver(X1[-nnn],X2[-nnn],X1[-1],X2[-1],c(0,theta),mesh,delt,nnn-1,T.seq[-nnn],strts,tro,secmom,seq1,seq2)
strts=tempp$like[,2:3]*recycle+strts*(1-recycle)
errors[i] = tempp$max
lnew=sum(tempp$like[,1])
rat=min(exp(lnew-lold)*pp(theta)/pp(theta.temp),1)
if(is.na(rat)){retry.count=retry.count+1}
if(retry.count>10){print('Fail: Local retry fail!');failed.chain=T;break;break;}
}
decodes=dec.matrix[,i]
}
u=runif(1)
is.true =(rat>u)
is.false=!is.true
theta=theta*is.true+theta.temp*is.false
lold=lnew*is.true +lold*is.false
par.matrix[,i]=theta
ll[i]=lold
probs[i] = mean(rs$probs)
kk=kk+is.true
acc[i]=is.true
if(decode)
{
ltemp1 = (tempp$like3+log(pmax(1-tempp$probs,0.000001)))*is.true +ltemp1*is.false
ltemp2 = (tempp$like2+log(pmax(tempp$probs,0.000001) ))*is.true +ltemp2*is.false
dc.new = sample(c(0,1),nnn-1,replace=TRUE)
dlike.old = ltemp1*(decodes==1)+(ltemp2)*(decodes==0)
dlike.new = ltemp1*(dc.new==1) +(ltemp2)*(dc.new==0)
is.switch = (pmin(exp(dlike.new-dlike.old),1)>runif(nnn-1))
decodes = dc.new*is.switch+decodes*(!is.switch)
dec.matrix[,i] = decodes
if(any(is.na(decodes))){print('Fail: Decoding structure NAs.');failed.chain=TRUE;break;}
}
if(max.retries>5000){print('Fail: Failed evaluation limit exceeded!');failed.chain=T;break;}
if(any(is.na(theta))){print('Fail: Samples were NA! ');failed.chain=T;break;}
setTxtProgressBar(pb, i)
i=i+1
}
close(pb)
acc = cumsum(acc)/(1:updates)
}
if(adapt!=0)
{
pb <- txtProgressBar(1,updates,1,style = 1,width = 65)
failed.chain=F
for(i in 2:updates)
{
theta.temp=theta
if(i>min(5000,round(burns/2)))
{
theta=theta+(1-adapt)*rnorm(length(theta),sd=sqrt(2.38^2/length(theta)*diag(covvec)))+adapt*rnorm(length(theta),sd=sqrt(0.1^2/length(theta)*diag(covvec)))
}else
{
theta=theta+rnorm(length(theta),sd=sds)
}
prop.matrix[,i] = theta
tempp=solver(X1[-nnn],X2[-nnn],X1[-1],X2[-1],c(0,theta),mesh,delt,nnn-1,T.seq[-nnn],strts,tro,secmom)
strts=tempp$like[,2:3]
errors[i] = tempp$max
lnew=sum(tempp$like[-excl,1])
rat=min(exp(lnew-lold)*pp(theta)/pp(theta.temp),1)
u=runif(1)
is.true =(rat>u)
is.false=!is.true
theta=theta*is.true+theta.temp*is.false
lold=lnew*is.true +lold*is.false
par.matrix[,i]=theta
ll[i]=lold
kk=kk+is.true
acc[i]=kk/i
muvec=muvec +1/(i)*(theta-muvec)
covvec = covvec+1/(i)*((theta-muvec)%*%t(theta-muvec)-covvec)
if(any(is.na(theta))){print('Fail: NaN thetas observed.');failed.chain=T;plot.chain=F;break;}
setTxtProgressBar(pb, i)
}
close(pb)
}
tme.eval = function(start_time)
{
start_time = as.POSIXct(start_time)
dt = difftime(Sys.time(), start_time, units="secs")
format(.POSIXct(dt,tz="GMT"), "%H:%M:%S")
}
tme=tme.eval(tme)
theta =apply(par.matrix[,-c(1:burns)],1,mean)
meanD=mean(-2*ll[-c(1:burns)])
rs=solver(X1[-nnn],X2[-nnn],X1[-1],X2[-1],c(0,theta),mesh,delt,nnn-1,T.seq[-nnn],strts,tro,secmom,seq1,seq2)
pd=meanD-(-2*sum(rs$like[-excl,1]))
DIC=meanD+pd
actual.p=length(theta)
model.inf=list(elapsed.time=tme,time.homogeneous=c('Yes','No')[2-homo],p=actual.p,DIC=DIC,pd=pd,N=length(X[,1])-length(excl)+1,Tag=Tag,burns=burns)
Info2=c( buffer7,
paste0("Chain Updates : ",updates),
paste0("Burned Updates : ",burns),
paste0("Time Homogeneous : ",c('Yes','No')[2-homo]),
paste0("Data Resolution : ",c(paste0('Homogeneous: dt=',round(max(diff(T.seq)[-excl]),4)),paste0('Variable: min(dt)=',round(min(diff(T.seq)[-excl]),4),', max(dt)=',round(max(diff(T.seq)[-excl]),4)))[2-homo.res]),
paste0("# Removed Transits. : ",c("None",length(excl))[2-is.null(exclude)]),
paste0("Density approx. : ",sol.state),
paste0('Elapsed time : ',tme),
buffer3,
paste0("dim(theta) : ",round(actual.p,3)),
paste0("DIC : ",round(DIC,3)),
paste0("pd (eff. dim(theta)): ",round(pd,3)),
buffer1)
Info2=data.frame(matrix(Info2,length(Info2),1))
colnames(Info2)=''
if(print.output)
{
print(Info2,row.names = FALSE,right=F)
}
if(plot.chain)
{
nper=length(theta)
if(nper==1){par(mfrow=c(1,2))}
if(nper==2){par(mfrow=c(2,2))}
if(nper==3){par(mfrow=c(2,2))}
if(nper>3)
{
d1=1:((nper)+2)
d2=d1
O=outer(d1,d2)
test=O-((nper)+1)
test[test<0]=100
test=test[1:4,1:4]
test
wh=which(test==min(test))
d1=d1[col(test)[wh[1]]]
d2=d2[row(test)[wh[1]]]
par(mfrow=c(d1,d2))
}
if(palette=='mono')
{
cols =rep('#222299',nper)
}else{
cols=rainbow_hcl(nper, start = 10, end = 275,c=100,l=70)
}
ylabs=paste0('theta[',1:nper,']')
for(i in 1:nper)
{
plot(prop.matrix[i,],col='gray90',type='s',main=ylabs[i],xlab='Iteration',ylab='')
lines(par.matrix[i,],col=cols[i],type='s')
abline(v=burns,lty='dotdash')
if(adapt!=0){abline(v=min(5000,round(burns/2)),lty='dotted',col='red')}
}
plot(acc,type='l',ylim=c(0,1),col='darkblue',main='Accept. Rate',xlab='Iteration',ylab='%/100')
abline(h=seq(0,1,1/10),lty='dotted')
abline(v=burns,lty='dotdash')
abline(h=0.4,lty='solid',col='red',lwd=1.2)
abline(h=0.2,lty='solid',col='red',lwd=1.2)
lines(probs,col='purple',lty='dotted',lwd=2)
box()
if(decode)
{
plot(term.discont~time[-1],type='n',col='darkblue',main='Jump detection prob.',xlab='Time (t)',ylab='Prob.',ylim=c(0,1))
abline(h=seq(0,1,1/10),lty='dotted',col='gray75')
lines(c(rowSums(dec.matrix)/updates)~c(time[-1]-0.5*diff(time)[1]),col='steelblue',type='h')
}
}
ret=list(par.matrix=t(par.matrix),acceptance.rate=acc,elapsed.time=tme,model.info=model.inf,failed.chain=failed.chain,covvec=covvec,prop.matrix=t(prop.matrix),errors=errors,decode.prob=rowSums(dec.matrix[,burns:updates])/(updates-burns),dec.matrix=dec.matrix)
class(ret) = 'JGQD.mcmc'
return(ret)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.