############################################################################################### # This R file covers all the codes correspond to the figures and results # of Sections (4.1), (4.2), (4.3) and (4.4) of the paper # Maryam Zolghadr and Sergei Zuyev "Optimal design of dilution # experiments under volume constraints" available from # arXiv:1212.3151 [math.ST]. Note that equations' number are identical # to those in the paper. ############################################################################################### # Call the nessessary libraries # The liberary we use throughout the codes is R-library "medea". It can be downloaded from # the webpage of one the authurs Sergei Zuyev: http://www.math.chalmers.se/~zuev/download.html ## require(medea) ############################################################################################### # These are the main variables troughout the code and should be defined in advance. Note that # some of these variables might take new value several times. # n =: total number of available mice # x =: Phase space, we are looking for optimum doses on this range. We might make this range # narrow to draw some figures more clear. # spacing =: We discretized the phase space $x$, and spacing shows the step to generate the # sequence of $x$. # tol =: Parameter of go.renorm(medea) and optimize functions. # step =: Parameter of go.renorm(medea) function # lam or lambda=: Parameter of $\lambda$ appears almost in all equations in the paper. # u =: The upper bound of Uniform distribution. This mainly correponds to Equations # (10), (11) and (12) where we assume $\lambda\sim \Uniform(1,u)$. # a or alpha =: Shape parameter of Gamma distribution. This mainly correponds to Equations # (10), (11) and (12) where we assume $\lambda\sim \Gamma(\alpha,\beta)$. # beta =: Scale parameter of Gamma distribution where $\lambda\sim \Gamma(\alpha,\beta)$. # q =: The integral $int_0^{\infty}$ over $\lambda$ in Equations (10), (11) and (12), where # $\lambda\sim \Gamma(\alpha,1)$ cannot be solved. Thus, we cut Gamma distribution from # above with the error=0.001. # err =: Error we have committed in using truncated distribution instead of full distribution. # d =: Represents the coefficient of $G_1(\mu;\lambda)$ in Equation (57). # c1 =: Represents the coefficient of $T_1(\mu;\lambda)$ in Equation (57). # c21=c20=c2 =: Represents the coefficient of $T_2(\mu;\lambda)$ in Equation (57). # p =: Parameter in Equation (66). It can be between $0$ and $1$ as it is a probability. # lam1 =: Parameter in Equation (66) # lam2 =: Parameter in Equation (66) ## # Initial values of the parameters: n=30 spacing=0.001 x=seq(0.001,1,by=spacing) #Note that $0$ is excluded. tol=10^-10 step=10 lam=20 u=120 a=100 beta=1 err=0.001 q=qgamma(err,shape=a,scale=beta,lower.tail=FALSE,log.p=FALSE) d=1 c1=10^-4 c21=1 c20=1 p=0.05 lam1=25 lam2=150 ############################################################################################### # We use this function throughout the code alot. So we can keep it in the beginning. r=function(y) exp(-y)/(1-exp(-y))*y^2 # Define $r(y) in equation (26) $ ############################################################################################### # Section 4.1, Figure (1) # ############################################################################################### # Figure (1), Plot (a): Draws $r(y)$ in equation (26) ## y=seq(0,8,by=0.01) y0=1.59362 # $y0$ is $y_{max}$. This is the point where $r(y)$ attains its maximum ## # r=function(y) exp(-y)/(1-exp(-y))*y^2 # Define $r(y) in equation (26) $ ## plot(c(0,8),c(0,0.8),ty="n",xlab="",ylab="") lines(y,r(y),lwd=2) segments(y0,0,y0,r(y0),lty=2) ## dev.print(device=pdf,file="r-graph.pdf") ############################################################################################### # Figure (1), Plot (b): Draws the optimal dose volume for $n=30$ mice as a function of $\lambda$. ## z=seq(y0*n,200,len=100) # Generates sequence of $lambda for $x_max\leq 1/n$ ## plot(c(0,200),c(0,0.04),ty="n",xlab="",ylab="") segments(0.001,1/n,y0*n,1/n,lwd=2) lines(z,y0/z,lwd=2) abline(h=0) ## dev.print(device=pdf,file="dose-n30.pdf") ############################################################################################### # Figure (1), plot (c): Draws $g_1$ in equation (25) and the line satisfies theorem (2), in # order to show the optimum atoms correspond to the optimum measure when $\lambda=20$. ## x=seq(0.001,0.2,by=0.01) ## y=r(lam*x)/lam^2 # Gradient function $g_1$ in equation (25), when $lambda=20$ ## plot(c(0,0.2),c(0,0.002),ty="n",xlab="",ylab="") lines(x,y,lwd=2) ## x0=1/n y0=r(lam*x0)/lam^2 # $(x0,y0)$ is the point of maximum for $g_1(x;20)$ segments(x0,0,x0,y0,lty=2) abline(h=0) ## dr=function(x) exp(-x)/(1-exp(-x))*x*(2-x/(1-exp(-x))) # Draws a tangent line to the curve of ## # gradient function at the atom of optimum. t=seq(0,0.07,len=10) dy=y0+dr(x0*lam)/lam*(t-x0) lines(t,dy) points(x0,y0,pch=19) ## dev.print(device=pdf,file="tan-n30-lam20.pdf") ############################################################################################### # Figure (1), plot (d): Draws $g_1$ in equation (25) and the line satisfies theorem (2), in # order to show atoms when $\lambda=100$. ## lam=100 x=seq(0.001,4/lam,len=100) x0=1.59362/lam y0=r(lam*x0)/lam^2 # $(x0,y0)$ is the point of maximum for $g_1(x;100)$ ## y=r(lam*x)/lam^2 # Gradient function $g_1$ in equation (25), when $lambda=100$ ## plot(c(0,4/lam),c(0,8*10^(-5)),ty="n",xlab="",ylab="") lines(x,y,lwd=2) ## segments(x0,0,x0,y0,lty=2) ## segments(1/n,0,1/n,8*10^(-5),lty=3) ## abline(h=c(0,y0)) # Draws a tangent line to the curve of gradient function at the atom of optimum. points(x0,y0,pch=19) ## dev.print(device=pdf,file="tan-n30-lam100.pdf") ############################################################################################### # Section 4.2, Figure (2) # # Section 4.2, Uniform prior distribution, Corresponding results for the criterion function # # $G_2$. # ############################################################################################### # Draws $(u-1)*g_2(x;u)$ for $u=10,20,200,1000$. ## x=seq(0.001,1,by=spacing) ## tg2=function(x,u) x*(log(1-exp(-u*x))-log(1-exp(-x))) # Define $(u-1)*g_2(x;u)$ as a function of ## # $x$ and $u$, where $g_2(x;u)$ is the ## # equation (28). plot(c(0,1),c(0,0.6),ty="n",xlab="",ylab="") lines(x,tg2(x,10),lwd=2) lines(x,tg2(x,20),lwd=2) lines(x,tg2(x,200),lwd=2) lines(x,tg2(x,1000),lwd=2) ## dev.print(device=pdf,file="g2-unif.pdf") ############################################################################################### # The aim of this code is to find a threshold for $u$ named $u^*$ , where for $u\leq u^*$ # we need to take all the substrate, and a proportion of the substrate otherwise. ## G2=function(x,u) # Define $G_2(n\delta_{x};u)$ as a function of $x$ and $u$ where { # $G_2(\mu;u)$ is Equation (27). num=1-exp(-u*x) denum=1-exp(-x) r=n*x/(u-1)*log(num/denum) return(r) } G2=Vectorize(G2) ## xmax=function(u) # For a given $u$ optimises $G_2(n\delta_{x};u)$ over $x$. { # and returns the point where $G_2(n\delta_{x};u)$ attains ff=function(x) G2(x,u) # its maximum. return(optimise(ff,c(0,1),maximum=TRUE,tol=tol)$maximum) } xmax=Vectorize(xmax) ## # We tried several values of $u=2,10,100,1000$, as argument of xmax. It turns out that # regardless of the upper bound $u$ (small,big) $xmax(u)\geq 1/30$ for any given $u$. # Thus for any given $u$ we have to take all the substrate and divide it between $n=30$ mice. ############################################################################################### # Section 4.2, Figure (3) # # Section 4.2, Uniform prior distribution, Corresponding results for the criterion # # function $G_3$. # ############################################################################################### # The code applies R-library "medea" to find an optimal solution to a measure optimization # problem, when $G_3$ in Eq. (29) is a goal function, and under constraints (8) and (9). # Here we assume $\lambda\sim \Uniform(1,u)$. The result of this code shows that for any # given $u$, the optimal measure contains only one atom. ## # To use this code someone need to specify spacing, phase space $x$, parameter of prior # distribution $u$, number of mice $n$, and define goal function and its gradient function # as we did. Note that the constraints we have in our problem were included in the library # of medea. However, someone can define more constraints and then call them in go.renorm. # For more details read the help of go.renorm {medea}. ## # x=seq(0.001,1,by=spacing) # u=120 # Run the code for u=120 to get the LHS plots of Figure (3). # u=20 # Run the code for u=20 to get the RHS plots of Figure (3). tol=0.0001 ## G3.func=function(mu,x) # $G_3$ goal function. { flambda=function(lambda) { ex=exp(-lambda*x) temp=sum(x^2*ex/(1-ex)*mu) return(1/(u-1)*log(temp)) } res=integrate(Vectorize(flambda),1,u)$value return(res) } ## G3.grad=function(mu,x) # $G_3$ Gradient function. { flambda=function(lambda) { ex=exp(-lambda*x) return(sum(x^2*ex/(1-ex)*mu)) } ii=NULL for (xx in x) { fl=function(lambda) { exx=exp(-lambda*xx) return(xx^2*exx/(1-exx)/((u-1)*flambda(lambda))) } ii=c(ii,integrate(Vectorize(fl),1,u)$value) } return(ii) } ## # Run "go.renorm" to find an optimum measure and value of $G_3$ for any given $b$, # where $b=\int x\mu(dx)$. # As medea has not been programmed for inequality constraints yet, and the constraint on a # volume has a form of inequality, we need to run medea for any given $b$ in a range $(0,1]$. # Then find $b$ that maximizes $G_3$. ## GK=function(b) { res<-go.renorm(mu=NA,G3.func,G3.grad,minim=FALSE,tol=tol,space=x,No.steps=step,details=FALSE,visualize=FALSE,constr=bary,cval=c(n,b)) while(!res$opt) res<-go.renorm(res) # go.renorm is searching for an optimum measure, as long as the optimum attains. return(invisible(list(fval=res$fval,mes=res$mes,b=b,space=x))) } gk=function(b) GK(b)$fval m=optimize(gk,c(0,1),maximum=TRUE) b=m$maximum #An optimum $b$. ## # Run "go.renorm" for the optimum $b$ to find an optimum measure including atoms and masses # attached to them, and then plot the results. ## res<-go.renorm(mu=NA,G3.func,G3.grad,minim=FALSE,tol=tol,space=x,No.steps=step,details=FALSE,visualize=FALSE,constr=bary,cval=c(n,b)) while(!res$opt) res<-go.renorm(res) summary(res) #Shows the optimum atoms and masses attach to them. #Answer: mu(0.025)=29.9695 mu(0.026)=0.0305298 plot(res) dev.print(device=pdf,file="g3-medea-unif-u-120.pdf") # dev.print(device=pdf,file="g3-medea-unif-u-20.pdf") ############################################################################################### # Section 4.2, Figure (4) # # Section 4.2, Uniform prior distribution, Corresponding results for the criterion # # function $G_3$. # ############################################################################################### # The aim is to determine $u^*$, in Eq. (33). ## tol=10^(-10) ## tG3=function(x,u) # Defines $\tilde{G_3(x;u)}$ in the Equation right after Equation (32) { tG3=function(lam) log(r(lam*x)) return(integrate(tG3,1,u)$value/(u-1)) } tG3=Vectorize(tG3) ## xmax=function(u) #optimise $tG3$ over $x$. { ff=function(x) tG3(x,u) return(optimise(ff,c(0,1),maximum=TRUE,tol=tol)$maximum) } xmax=Vectorize(xmax) ## xm1=function(a) (xmax(a)-1/n)^2 optimize(xm1,c(50,100),tol=tol)$minimum # Gives solution to $x_max(u^*)=1/n$, so that $u^*=90.66$. ############################################################################################### # Figure (4) LHS: Draws $\tilde{G_3}$ in Equation right after Equation (32) for $u=50,100,200$. ## x=seq(0.001,0.2,by=0.001) tg200=tG3(x,200) tg100=tG3(x,100) tg50=tG3(x,50) ## plot(c(0,0.2),c(-15,0),ty="n",xlab="",ylab="") lines(x,tg50,lwd=2) lines(x,tg100,lwd=2) lines(x,tg200,lwd=2) abline(v=1/30,lty=2) ## dev.print(device=pdf,file="g3-unif-oneatom.pdf") ############################################################################################### # Figure (4) RHS: Draws $argmax(tilde{G_3})$ for $n=30$ mice as a function of $u$. xm=xmax(10:200) plot(c(10,200),c(0,0.3),ty="n",xlab="u",ylab="x_max") lines(10:200,xm,lwd=2) abline(h=1/30,lty=2) abline(h=0) abline(h=seq(0.01,0.3,by=0.01),col="lightgray",lty=3) abline(v=seq(10,200,by=10),col="lightgray",lty=3) dev.print(device=pdf,file="g3-unif-xmax.pdf") ############################################################################################### # Section 4.2, Uniform prior distribution, Corresponding results for the criterion function # # $G_4$. # # Section 4.2, Figure (5). # ############################################################################################### # The code applies R-library "medea" to find an optimal solution to a measure optimization # problem, when $G_4$ is the Eq. (34) and a goal function, and under constraints (8) and (9). # Here we assume $\lambda\sim \Uniform(1,u)$. The result of this code shows that for any # given $u$, the optimal measure contains only one atom. Run it for $u=100$ and $u=20$. ## # u=120 # u=20 x=seq(0.001,1,by=spacing) tol=0.0001 ## G4.func=function(mu,x) # $G_4$ goal function. { flambda=function(lambda) { ex=exp(-lambda*x) temp=sum(x^2*ex/(1-ex)*mu) return(1/((u-1)*temp)) } res=integrate(Vectorize(flambda),1,u)$value return(-res) } ## G4.grad=function(mu,x) # $G_4$ Gradient function. { flambda=function(lambda) { ex=exp(-lambda*x) temp=sum(x^2*ex/(1-ex)*mu) return(temp^2) } ii=NULL for (xx in x) { fl=function(lambda) { exx=exp(-lambda*xx) return(xx^2*exx/(1-exx)/((u-1)*flambda(lambda))) } ii=c(ii,integrate(Vectorize(fl),1,u)$value) } return(ii) } ## G4.GK=function(b) # Find an optimum $b$. { res<-go.renorm(mu=NA,G4.func,G4.grad,minim=FALSE,tol=tol,space=x,No.steps=step,details=FALSE,visualize=FALSE,constr=bary,cval=c(n,b)) while(!res$opt) res<-go.renorm(res) return(invisible(list(fval=res$fval,mes=res$mes,b=b,space=x))) } G4.gk=function(b) G4.GK(b)$fval G4.m=optimize(G4.gk,c(0,1),maximum=TRUE) G4.b=G4.m$maximum ## # Use the optimum $b$ to find an optimum measure. ## G4.res<-go.renorm(mu=NA,G4.func,G4.grad,minim=FALSE,tol=tol,space=x,No.steps=step,details=FALSE,visualize=FALSE,constr=bary,cval=c(n,G4.b)) while(!G4.res$opt) res<-go.renorm(G4.res) summary(G4.res) #mu(0.018)=29.9866 mu(0.019)=0.0134 plot(G4.res) dev.print(device=pdf,file="G4-medea-unif-u-120.pdf") # dev.print(device=pdf,file="G4-medea-unif-u-20.pdf") ############################################################################################### # The aim of this part is to specify $u^*$, in Eq. (33), where the optimization # problem is maximisation of $G_4(n*\delta_{x};u)$ such that $x^*\leq 1/30$, where $x^*$ is # the representative of the optimum dose. ## tol=10^(-10) ## G4=function(x,u) # Define $n*G_4(n\delta_{x};u)$ in equation (36) { num=exp(u*x)-exp(x)-x*(u-1) den=x^3*(u-1) return(-num/den) } G4=Vectorize(G4) ## xmax=function(u) #optimise $G_4$ over $x$. { ff=function(x) G4(x,u) return(optimise(ff,c(0,1),maximum=TRUE,tol=tol)$maximum) } xmax=Vectorize(xmax) ## xm=function(a) (xmax(a)-1/n)^2 optimize(xm,c(20,100),tol=tol)$minimum #Gives solution to x_max(u^*)=1/30 so that u^*= 64.47. ############################################################################################### # Figure (5) LHS ## x=seq(0.001,0.2,by=0.001) g100=-log(abs(G4(x,100))) g50=-log(abs(G4(x,50))) g20=-log(abs(G4(x,20))) plot(c(0.0,0.2),c(-12,-4),ty="n",xlab="",ylab="") lines(x,g20,lwd=2) lines(x,g50,lwd=2) lines(x,g100,lwd=2) abline(v=1/30,lty=2) dev.print(device=pdf,file="g4-unif-oneatom.pdf") ############################################################################################### # Figure (5) RHS ## xmm=xmax(10:200) plot(c(10,200),c(0,0.2),ty="n",xlab="u",ylab="x_max") lines(10:200,xmm,lwd=2) abline(h=1/n,lty=2) abline(h=0) abline(h=seq(0.02,0.2,by=0.02),col="lightgray",lty=3) abline(v=seq(10,200,by=10),col="lightgray",lty=3) dev.print(device=pdf,file="g4-xmax.pdf") ############################################################################################### # Section 4.2, Figure (6). # # Section 4.2, Gamma prior distribution, Corresponding results for the criterion function # # $G_2$. # ############################################################################################### # The aim of this part is to specify $u^*$, in Eq. (42) in our paper, where the optimization # problem is maximisation of $G_2(n\delta_{x};\alpha,\beta)$ in Equation (39) such that # $x^* \leq 1/n$. ## g2=function(x,alpha,beta=1,err=10^(-3)) # Defines gradient function $g_2(x,\mu;\alpha,\beta)$ { # in Equation (41). f=function(lam) r(lam)*dgamma(lam,alpha-2,rate=beta/x) return(integrate(f,qgamma(err,alpha-2,rate=beta/x),qgamma(1-err,alpha-2,rate=beta/x),abs.tol=err)$value/(alpha-1)/(alpha-2)*beta^2) } g2=Vectorize(g2) ## xmax=function(alpha,beta) # optimise $g_2$ over $x$. { ff=function(x) g2(x,alpha,beta=1) return(optimise(ff,c(0,1),maximum=TRUE,tol=tol)$maximum) } xmax=Vectorize(xmax) ## xm1=function(a) xmax(a,1) xm1=function(a) (xmax(a,1)-1/n)^2 optimize(xm1,c(20,70),tol=tol)$minimum # Gives solution to $x_max(\alpha^*)=1/n$ so that $\alpha^*= 49.68$. ############################################################################################### # Figure (6) LHS. ## plot(c(0,0.2),c(0,8*10^(-4)),ty="n",xlab="x",ylab="") abline(h=0) lines(x,g2(x,30),lwd=2) lines(x,g2(x,50),lwd=2) lines(x,g2(x,75),lwd=2) lines(x,g2(x,100),lwd=2) segments(1/n,0,1/n,0.0008,lty=2) ## dev.print(device=pdf,file="g2-gamma.pdf") ############################################################################################### # Figure (6) RHS ## xm=xmax(10:200,beta=1) plot(c(10,200),c(0,0.2),ty="n",xlab="alpha",ylab="x_max") lines(10:200,xm,lwd=2) abline(h=1/n,lty=2) abline(h=0) abline(h=seq(0.01,0.2,by=0.01),col="lightgray",lty=3) abline(v=seq(10,200,by=10),col="lightgray",lty=3) ## dev.print(device=pdf,file="g2-xmax-gamma.pdf") ############################################################################################### # Section 4.2, Gamma prior distribution, Corresponding results for the criterion function # # $G_3$. # # Section 4.2, Figure (7). # ############################################################################################### # Applies R-library "medea" to find an optimal solution to a measure optimization problem, # when $G_3$ in Eq. (43) is a goal function, and under constraints (8) and (9). Here we assume # $\lambda\sim \Gamma(\alpha,1)$. The result of this code shows that for any given $\alpha$, # the optimal measure contains only one atom. ## # a=100 # a=20 x=seq(0.001,1,by=spacing) tol=0.0001 ## G3.func=function(mu,x) # $G_3$ goal function { flambda=function(lambda) { ex=exp(-lambda*x) temp=sum(x^2*ex/(1-ex)*mu) r1=log(temp) r2=dgamma(lambda,shape=a,scale=beta,log=FALSE) r3=r1*r2 return(r3) } res=integrate(Vectorize(flambda),0,q)$value return(res) } ## G3.grad=function(mu,x) # $G_3$ Gradient function { flambda=function(lambda) #denumerator { ex=exp(-lambda*x) return(sum(x^2*ex/(1-ex)*mu)) } ii=NULL for (xx in x) { fl=function(lambda) { exx=exp(-lambda*xx) r1=(xx^2*exx)/((1-exx)*flambda(lambda)) r2=dgamma(lambda,shape=a,scale=beta,log=FALSE) r3=r1*r2 return(r3) } ii=c(ii,integrate(Vectorize(fl),0,q)$value) } return(ii) } ## GK=function(b) # Find an optimum $b$. { res<-go.renorm(mu=NA,G3.func,G3.grad,minim=FALSE,tol=tol,space=x,No.steps=step,details=FALSE,visualize=FALSE,constr=bary,cval=c(n,b)) while(!res$opt) res<-go.renorm(res) return(invisible(list(fval=res$fval,mes=res$mes,b=b,space=x))) } gk=function(b) GK(b)$fval m=optimize(gk,c(0,1),maximum=TRUE) b=m$maximum ## # Use the optimum $b$ to find an optimum measure. ## res<-go.renorm(mu=NA,G3.func,G3.grad,minim=FALSE,tol=tol,space=x,No.steps=step,details=FALSE,visualize=FALSE,constr=bary,cval=c(n,b)) while(!res$opt) res<-go.renorm(res) summary(res) #mu(0.015)=0.0227 mu(0.016)=29.9773 plot(res) ## dev.print(device=pdf,file="G3-medea-Gamma-a-100.pdf") #dev.print(device=pdf,file="G3-medea-Gamma-a-20.pdf") ############################################################################################### # The aim of this part is to specify $a^*$, in Eq. (42) in our paper, where the optimization # problem is maximisation of $G_3(n\delta_{x};\alpha,\beta)$ in Equation (43) such that # $x^* \leq 1/n$. ## tol=10^-10 ## # Defines $tG3=G_3(n\delta_{x};\alpha,\beta)-c(\alpha,\beta)$, where # $G_3(n\delta_{x};\alpha,\beta)$ is the Equation (43) for $\mu=n\delta_{x}$, and # $c(\alpha,\beta)=\int_{0}^{\infty}f_{\alpha,\beta}(\lambda)\logn$ does not depend on $x$. # Note that $n=30$ mice. tG3=function(x,alpha,beta=1,err=10^(-3)) { f=function(lambda) log(r(lambda*x))*dgamma(lambda,shape=alpha,scale=beta,log=FALSE) return(integrate(f,0,qgamma(err,shape=alpha,rate=beta,lower.tail=FALSE,log.p=FALSE))$value) } tG3=Vectorize(tG3) ## xmax=function(a) ##optimise $G_3$ over $x$. { ff=function(x) tG3(x,a) return(optimise(ff,c(0,1),maximum=TRUE,tol=tol)$maximum) } xmax=Vectorize(xmax) ## xm=function(a) (xmax(a)-1/n)^2 optimize(xm,c(20,100),tol=tol)$minimum #Gives solution to $x_max(a^*)=1/n$ so that $a^*=47.70 ############################################################################################### #Figure (7) LHS ## x=seq(0.001,0.5,by=0.001) g200=tG3(x,200) g100=tG3(x,100) g50=tG3(x,50) plot(c(0,0.5),c(-14,0),ty="n",xlab="",ylab="") #-15 lines(x,g50,lwd=2) lines(x,g100,lwd=2) lines(x,g200,lwd=2) abline(v=1/n,lty=2) ## dev.print(device=pdf,file="g3-Gamma-oneatom.pdf") ############################################################################################### #Figure (7) RHS ## xm1=xmax(10:200) plot(c(10,200),c(0, 0.16),ty="n",xlab="alpha",ylab="x_max") lines(10:200,xm1,lwd=2) abline(h=1/n,lty=2) abline(h=0) abline(h=seq(0.01,0.16,by=0.01),col="lightgray",lty=3) abline(v=seq(10,200,by=10),col="lightgray",lty=3) ## dev.print(device=pdf,file="g3-Gamma-xmax.pdf") ############################################################################################### # Section 4.2, Gamma prior distribution, Corresponding results for the criterion function # $G_4$. # Section 4.2, Figure (8). ############################################################################################### # Applies R-library "medea" to find an optimal solution to a measure optimization problem, # when $G_4$ in Eq. (45) is a goal function, and under constraints (8) and (9). Here we # assume $\lambda\sim \Gamma(\alpha,1)$. The result of this code shows that for any given # $\alpha$, the optimal measure contains only one atom. ## x=seq(0.001,1,by=spacing) tol=0.0001 ## G4.func=function(mu,x) # $G_4$ goal function { flambda=function(lambda) { ex=exp(-lambda*x) r1=sum(x^2*ex/(1-ex)*mu) r2=dgamma(lambda,shape=a,scale=beta,log=FALSE) r3=r2/r1 return(r3) } res=integrate(Vectorize(flambda),0,q)$value return(-res) } ## G4.grad=function(mu,x) # $G_4$ Gradient function { flambda=function(lambda) { ex=exp(-lambda*x) temp=sum(x^2*ex/(1-ex)*mu) return(temp^2) } ii=NULL for (xx in x) { fl=function(lambda) { exx=exp(-lambda*xx) r1=xx^2*exx/((1-exx)*flambda(lambda)) r2=dgamma(lambda,shape=a,scale=beta,log=FALSE) r3=r1*r2 return(r3) } ii=c(ii,integrate(Vectorize(fl),0,q)$value) } return(ii) } ## G4.GK=function(b) # Find an optimum $b$ { res<-go.renorm(mu=NA,G4.func,G4.grad,minim=FALSE,tol=tol,space=x,No.steps=step,details=FALSE,visualize=FALSE,constr=bary,cval=c(n,b)) while(!res$opt) res<-go.renorm(res) return(invisible(list(fval=res$fval,mes=res$mes,b=b,space=x))) } G4.gk=function(b) G4.GK(b)$fval G4.m=optimize(G4.gk,c(0,1),maximum=TRUE) G4.b=G4.m$maximum ## # Use the optimum $b$ to find an optimum measure. ## G4.res<-go.renorm(mu=NA,G4.func,G4.grad,minim=FALSE,tol=tol,space=x,No.steps=step,details=FALSE,visualize=FALSE,constr=bary,cval=c(n,G4.b)) while(!G4.res$opt) res<-go.renorm(G4.res) summary(G4.res) #mu(0.015)=0.0309069 mu(0.016)=29.9691 plot(G4.res) ## dev.print(device=pdf,file="G4-medea-Gamma-a-100.pdf") #dev.print(device=pdf,file="G4-medea-Gamma-a-20.pdf") ############################################################################################### # The aim of this part is to specify $a^*$, in Eq. (42) in our paper, where the optimization # problem is maximisation of $G_4$ in Equation (47) such that $x^* \leq 1/n$. ## tol=10^-10 ## G4=function(x,alpha,beta=1,err=10^(-3)) # The code defines $n G_4(n\delta_{x};\alpha)$ in Equation (47). { f=function(lambda) -(1-exp(-lambda*x))/(x^2*exp(-lambda*x))*dgamma(lambda,shape=alpha,scale=beta,log=FALSE) return(integrate(f,0,qgamma(err,shape=alpha,rate=beta,lower.tail=FALSE,log.p=FALSE))$value) } G4=Vectorize(G4) ## xmax=function(a) #optimise $G_4$ over $x$. { ff=function(x) G4(x,a) return(optimise(ff,c(0,1),maximum=TRUE,tol=tol)$maximum) } xmax=Vectorize(xmax) ## xm=function(a) (xmax(a)-1/n)^2 optimize(xm,c(20,100),tol=tol)$minimum #Gives solution to x_max(a^*)=1/n so that a^*=45.74 ############################################################################################### # Figure (8) LHS ## x=seq(0.001,0.2,by=0.001) g100=-log(abs(G4(x,100))) g50=-log(abs(G4(x,50))) g20=-log(abs(G4(x,20))) ## plot(c(0.0,0.2),c(-18,-6),ty="n",xlab="",ylab="") lines(x,g20,lwd=2) lines(x,g50,lwd=2) lines(x,g100,lwd=2) abline(v=1/n,lty=2) dev.print(device=pdf,file="g4-Gamma-oneatom.pdf") ############################################################################################### # Figure (8) RHS ## xm1=xmax(10:200) plot(c(10,200),c(0,0.15),ty="n",xlab="alpha",ylab="x_max") lines(10:200,xm1,lwd=2) abline(h=1/n,lty=2) abline(h=0) abline(h=seq(0.02,0.2,by=0.02),col="lightgray",lty=3) abline(v=seq(10,200,by=10),col="lightgray",lty=3) ## dev.print(device=pdf,file="g4-Gamma-xmax.pdf") ############################################################################################### # Section 4.3, Figure (9). # # Section 4.3, Figure (10). # ############################################################################################### # Find an optimal dose when $\tilde{G_1}$ in Equation (57) is a goal function and an optimum # measure contains one atom. ## x=seq(0.001,0.05,by=spacing) ## G1.func=function(x,lambda) # $G1(\mu;\lambda)$ in Equation (7) where $\mu=n\delta_{x}$ { ex=exp(-lambda*x) G1=n*(x^2)*(ex/(1-ex)) return(G1) } G1.func=Vectorize(G1.func) ## T1.func=function(x,lambda) # $T_1(\mu;\lambda)$ in Equation (49) where $\mu=n\delta_{x}$ { ex=exp(-lambda*x) T1=ex*n return(T1) } ## T21.func=function(x,lambda) # T_{21}(\mu;\lambda) in Equation (52) where $\mu=n\delta_{x}$ { b=n*x T21=exp(-lambda*b) return(T21) } ## T20.func=function(x,lambda) # T_{20}(\mu;\lambda) in Equation (53) where $\mu=n\delta_{x}$ { ex=exp(-lambda*x) expon=n*log(1-ex) res=exp(expon) return(res) } T20.func=Vectorize(T20.func) ## tG1.func=function(x,lambda) # $\tilde{G_1}$ in Equation (57) where $\mu=n\delta_{x}$ { tG1=d*G1.func(x,lambda)-c1*T1.func(x,lambda)-c21*T21.func(x,lambda)-c20*T20.func(x,lambda) return(tG1) } ## xmax=function(lambda) # Optimize $\tilde{G_1}(n\delta_{x};\lambda)$ over $x$ { ff=function(x) tG1.func(x,lambda) return(optimise(ff,c(0,0.04),maximum=TRUE,tol=tol)$maximum) } xmax=Vectorize(xmax) ## xm1=function(lambda) (xmax(lambda)-1/n)^2 opt.xm1=optimize(xm1,c(20,100),tol=tol)$minimum #41.80 ############################################################################################### # Figure (9) ## x=seq(0.001,0.04,by=spacing) ## # Draws the solid line, the optimal dose volume for $n=30$ mice as a function of $\lambda$, # when $G_1(n\delta_{x};\lambda)$ is a goal function. This solid line is identical to the # plot (b) in Figure (1). x0=1.59362 z=seq(x0*n,200,len=100) plot(c(0,200),c(0,0.04),ty="n",xlab="",ylab="") segments(0,1/n,x0*n,1/n,lwd=2) lines(z,x0/z,lwd=2) abline(h=0) ## # Draws the dashed line, the optimal dose volume for $n=30$ mice as a function of $\lambda$, # when $\tilde{G_1}(\mu;\lambda)$ is a goal function and $\mu=n\delta_{x}$. ## lam=seq(opt.xm1,200,len=100) #values of lambda for which argmax(G1(n\delta_{x};lambda))<1/n xx=xmax(lam) lines(lam,xx,lty=2) abline(h=0) dev.print(device=pdf,file="tildeG1-does-lam.pdf") ############################################################################################### # Figure (10): Find the optimal dose when $\tilde{G_4}$ in Equation (64) is a goal function, # and $\mu=n\delta_{x}$ the optimum measure contains one atom. ## x=seq(0.001,0.8,by=spacing) c1=5 c21=10^4 c20=10^4 a=100 q=qgamma(err,shape=a,scale=beta,lower.tail=FALSE,log.p=FALSE) #c1=0.005 #c21=5 #c20=5 ## G4.func=function(x,alpha) # $G_4$ in Equation (47) { flambda=function(lambda) { gamma=dgamma(lambda,shape=alpha,scale=beta,log=FALSE) G4=-(1-exp(-lambda*x))/(n*x^2*exp(-lambda*x))*gamma return(G4) } return(integrate(flambda,0,q)$value) } G4.func=Vectorize(G4.func) ## T1.func=function(x,alpha) # $T_1$ in Equation (59) where $\mu=n\delta_{x}$ { T1=n/(1+x)^alpha return(T1) } ## T21.func=function(x,alpha) # $T_21$: the first term in the RHS of Equation (62) where { # $\mu=n\delta_{x}$. b=n*x T21=(1/(1+b))^alpha return(T21) } ## T20.func=function(x,alpha) # $T_21$: the second term in the RHS of Equation (62) where { # $\mu=n\delta_{x}$. flambda=function(lambda) { gamma=dgamma(lambda,shape=alpha,scale=beta,log=FALSE) ex=exp(-lambda*x) expon=n*log(1-ex) res=exp(expon)*gamma return(res) } return(integrate(flambda,0,q)$value) } T20.func=Vectorize(T20.func) ## tG4.func=function(x,alpha) # $tilde{G_4}$ in Equation (64) where $\mu=n\delta_{x}$. { tG4=d*G4.func(x,alpha)-c1*T1.func(x,alpha)-c21*T21.func(x,alpha)-c20*T20.func(x,alpha) return(tG4) } tG4.func=Vectorize(tG4.func) ## plot(c(0,0.08),c(-10,-4),ty="n",xlab="x",ylab="") abline(v=1/n,lty=2) alpha=30 #0.043 #0.05 g30=-log(abs(tG4.func(x,alpha))) lines(x,g30,lwd=2) ## alpha=50 #0.028 #0.031 g50=-log(abs(tG4.func(x,alpha))) lines(x,g50,lwd=2) ## alpha=100 #0.015 #0.016 g100=-log(abs(tG4.func(x,alpha))) lines(x,g100,lwd=2) dev.print(device=pdf,file="tildeG4-one-atom-Gamma.pdf") ############################################################################################### # Section 4.4, Figure (11). # ############################################################################################### # Applies R-library "medea" to find an optimal solution to a measure optimization problem, # when $G_{1M}$ in Eq. (66) is a goal function, and under constraints (8) and (9). For the # specific values of parameters $p_1=0.05$, $\lambda_1=25$ and $\lambda_2=150$ in Equation # (66), it turns out that the optimum measure contains 2 atoms. x=seq(0.01,0.1,by=spacing) # or x=seq(0.001,1,by=spacing) tol=0.0001 ## G1.func=function(mu,x,lambda) # G1 goal function in Equation (7) { ex=exp(-lambda*x) G1=sum(x^2*ex/(1-ex)*mu) return(G1) } ## G1.grad=function(mu,x,lambda) # G1 Gradient function in Equation (25) { ex=exp(-lambda*(x)) g1=x^2*(ex/(1-ex)) return(g1) } G1.grad=Vectorize(G1.grad) ## G1.M.func=function(mu,x,lambda) # G1.M goal function { G1.M=p*G1.func(mu,x,lam1)+(1-p)*G1.func(mu,x,lam2) return(G1.M) } ## G1.M.grad=function(mu,x,lambda) # G1.M Gradient function { G1.M=p*G1.grad(mu,x,lam1)+(1-p)*G1.grad(mu,x,lam2) return(G1.M) } ## G1.GK=function(b) # Find ab optimum $b$ { res<-go.renorm(mu=NA,G1.M.func,G1.M.grad,minim=FALSE,tol=tol,space=x,No.steps=step,details=FALSE,visualize=TRUE,constr=bary,cval=c(n,b)) while(!res$opt) res<-go.renorm(res) return(invisible(list(fval=res$fval,mes=res$mes,b=b,space=x))) } G1.gk=function(b) G1.GK(b)$fval G1.m=optimize(G1.gk,c(0,1),maximum=TRUE) G1.b=G1.m$maximum ## # Use the optimum $b$ to find an optimum measure. ## G1.res<-go.renorm(mu=NA,G1.M.func,G1.M.grad,minim=FALSE,tol=tol,space=x,No.steps=step,details=FALSE,visualize=TRUE,constr=bary,cval=c(n,G1.b)) while(!G1.res$opt) res<-go.renorm(G1.res) summary(G1.res) plot(G1.res) dev.print(device=pdf,file="g1M-medea.pdf") ############################################################################################### # The End # ###############################################################################################