#source("c:/Work07/ShortCourse07/Programs/PKPred.txt",print.eval=TRUE)#
#Does prediction for PK model - Individual
data<-read.table("c:/WorkHome06/ShortCourse06/Programs/PKPredData.txt",sep = "",header=T)
Sdata<-1/data[,1] 		#Winbugs give precision
Cdata<-data[,2:4]
Vdata<-data[,5:7]
t=c(2,4,6,8,10,24)			#times
g<-function(C,V,t){(30/V)*exp(-C*t/V)}
mC<-mean(Cdata)
mV<-mean(Vdata)
#-----------------Means for individual 1------------------------------------
m1<-g(mC[1],mV[1],2);m2<-g(mC[1],mV[1],6);m3<-g(mC[1],mV[1],10);s<-sqrt(mean(Sdata))
#-----------------Predictive Density------------------------------------------
pred<-function(x,t,C,V,S)
{
m<-length(C);temp<-0
for(i in 1:m)temp<-temp+dnorm(x,mean=g(C[i],V[i],t),sd=sqrt(S[i]))
return(temp/m)
}
#------------------Plot-----------------------------------------------
par(mfrow=c(1,3));ymax=6;d<-5
plot(function(x)dnorm(x, mean=m1,sd=s),xlim=c(m1-d*s,m1+d*s),ylim=c(0,ymax),main="Time=2",ylab="")
par(new=T)
plot(function(x)pred(x,2,Cdata[,1],Vdata[,1],Sdata),xlim=c(m1-d*s,m1+d*s),ylim=c(0,ymax),lty=2,ylab="")
plot(function(x)dnorm(x, mean=m2,sd=s),xlim=c(m2-d*s,m2+d*s),ylim=c(0,ymax),main="Time=6",ylab="")
par(new=T)
plot(function(x)pred(x,6,Cdata[,1],Vdata[,1],Sdata),xlim=c(m2-d*s,m2+d*s),ylim=c(0,ymax),lty=2,ylab="")
plot(function(x)dnorm(x, mean=m3,sd=s),xlim=c(m3-d*s,m3+d*s),ylim=c(0,ymax),main="Time=10",ylab="")
par(new=T)
plot(function(x)pred(x,10,Cdata[,1],Vdata[,1],Sdata),xlim=c(m3-d*s,m3+d*s),ylim=c(0,ymax),lty=2,ylab="")