Differentiability and approximation of probability functions under Gaussian mixture Models.

#this chunk is to install and require the knitr package, that will produce an html or pdf output report 
if (!require("knitr")) {
  install.packages("knitr") 
  require("knitr") 
}

## Loading required package: knitr

The setup of this example is that described in section 5 of the article “Differentiability and approximation of probability functions under Gaussian mixture Models” using parameters, sample sizes and other numerical specifications chosen arbitrarily. You may download the file GMMpf.Rmd and compile it using R to change these parameters to your taste.

The example

Suppose \[\xi \sim \mathcal{N}\left(\begin{pmatrix} 0 \\ C \end{pmatrix},\begin{bmatrix} 1-C^2 & 0 \\ 0& 1-C^2 \end{bmatrix}\right), \] where the y coordinate \(C \in [-1,1]\) is such that \(\frac{C+1}{2}\) has, for some \(\delta\in (0,1)\), the Beta(\(\delta,\delta\)) distribution, i.e. the density of random variable \(C\) is \(\frac{\Gamma(2\delta)}{2^{2\delta-1}\Gamma(\delta)^2}(1-c^2)^{\delta-1}\) for \(c\in [-1,1]\). The chance constraint \(g:\mathbb{R} \times \mathbb{R}^2\) is given by \[g(x,\xi)=||\xi||^2-x^2-2.\] It is clear that \(g(x,(0,c)^\top) \leq 0\) for each \(c\in [-1,1]\) and that, for \(r>0\) and \(v=(v_1,v_2)^\top\) such that \(v_1^2+v_2^2=1\), using the spherical-radial decomposition $=(0,c)^+r(c) v $ with \(r>0\) and \(\mathcal{L}(c)= \begin{bmatrix} \sqrt{1-C^2} & 0 \\ 0& \sqrt{1-C^2} \end{bmatrix}\) it holds that

\[ g(x,(0,c)^\top +r\mathcal{L}(c) v ) \leq 0 \iff r\leq \rho_c(x,v)=\frac{-cv_2+\sqrt{x^2+2-c^2v_1^2}}{\sqrt{1-c^2}}.\]

The radial probability function and its gradient are therefore given by \[e_c(x,v)= \int_0^{\rho_c(x,v)}re^{-r^2/2}dt=1-e^{-\rho_c(x,v)^2/2},\]

and \[\nabla_x e_c(x,v)=\frac{ x\rho_c(x,v)e^{-\rho_c(x,v)^2/2}}{ \sqrt{ [1-c^2] [ x ^2+2 - c^2v_1 ^2] } }\]

The probability function and its gradient are given by

\[\Phi(x) =\frac{\Gamma(2\delta)}{2^{2\delta-1}\Gamma(\delta)^2}\int_{-1}^1(1-c^2)^{\delta-1} F_\chi^{2,c^2/(1-c^2)}([1-c^2][x^2+2])dc,\] and \[\nabla \Phi (x)=\frac{x\Gamma(2\delta)}{2^{2\delta}\Gamma(\delta)^2}\int_{-1}^1(1-c^2)^{\delta-1} f_\chi^{2,c^2/(1-c^2)}([1-c^2][x^2+2])dc,\]

where \(F_\chi^{k,\lambda}\) denotes the cumulative distribution function of a non-central chi-squared random variable with \(k\) degrees of freedom and noncentrality parameter \(\lambda\), and \(f_\chi^{k,\lambda}\) denoting its corresponding density. We now conduct a simulation based estimation of \(\Phi\) and its gradient via

\[\Phi_N(x)=\frac{1}{N}\sum_{i=1}^N e_{c_i}(x,v_i) \quad \text{ and } \quad \nabla\Phi_N(x)=\frac{1}{N}\sum_{i=1}^N \nabla_x e_{c_i}(x,v_i).\] comparing \(\Phi_N\) with the estimator based on the proportion of simulations verifying the constraint

\[\tilde{\Phi}_N(x)=\frac{1}{N}\sum_{i=1}^N 1_{ \{ g(x,\xi_i)\leq 0 \}}.\]

Parameter setup

set.seed(681986) #seeded for replicability, change or remove for alternate samples
N<-100 #sample size
d<-0.5 #delta parameter, see more details at end
xlow<-0 #lower bound for range of decision variable x
xup<-4 #upper bound for range of decision variable x
xvalues<-200 #number of different values of x in range to simulate the probability function and estimators

Simulation of probability function curves and functional estimates

estradial<-rep(0,xvalues)
estgradient<-rep(0,xvalues)
estprop<-rep(0,xvalues)
trueprob<-rep(0,xvalues)
truegradient<-rep(0,xvalues)
v<-runif(N,0,2*pi)
c<-2*rbeta(N,d,d)-1
x<-seq(xlow,xup, length.out=xvalues)
rs<-rep(0,N)
for(w in 1:N){
  rs[w]<-rnorm(1,0,sqrt(1-c[w]^2))^2+rnorm(1,c[w],sqrt(1-c[w]^2))^2
}
for(i in 1:xvalues){
  rcrit<-(-c*sin(v)+sqrt(x[i]^2+2 -c^2*cos(v)^2))/sqrt(1-c^2)
  estradial[i]<-mean(1-exp(-rcrit^2/2))
  estgradient[i]<-mean(x[i]*rcrit*exp(-rcrit^2/2)/sqrt((x[i]^2+2 -c^2*cos(v)^2)*(1-c^2)))
  estprop[i]<-sum(rs < 2+x[i]^2)/N
  cdf<-function(y){(1-y^2)^(d-1)*pchisq((x[i]^2+2)/(1-y^2),2,ncp = y^2/(1-y^2))}
  trueprob[i]<-integrate(cdf,-1,1,abs.tol = 0.000001)$value
  dens<-function(y){(1-y^2)^(d-2)*dchisq((x[i]^2+2)/(1-y^2),2,ncp = y^2/(1-y^2))}
  truegradient[i]<-2*integrate(dens,0,1,abs.tol = 0.000000001)$value*x[i]
}

trueprob<-trueprob/(2^(2*d-1)*beta(d,d))
truegradient<-truegradient/(2^(2*d-2)*beta(d,d))

Plots

par(mfrow=c(1,2))
plot(x,trueprob,type="l", ylab = "Probability", main=paste("Constraint probability. N=", N))
lines(x,estprop,col="blue")
lines(x,estradial,col="red")
legend("bottomright", col = c("black", "red", "blue"), lty = 1, legend = c(expression(Phi * " (true)"), expression(Phi[N]* " (prop)"), expression(tilde(Phi)[N]* " (naive)")))
plot(x,truegradient,type="l", ylab = "Derivative", main="Gradient of probability function")
lines(x,estgradient,col="red")
legend("topright", col = c("black", "red"), lty = 1, legend = c(expression(paste(nabla,Phi* " (true)")), expression(paste(nabla,Phi[N]* " (prop)"))))

Visualization of mixing distributions

This part contains a visualization of the density of the vector \(\xi\) for different Beta parameter \(\delta\), both as a heatplot and as a 3 dimensional plot. Default is the value in the example \(\delta=0.5\). A value \(\delta \geq 1\) will not produce a bimodal distribution.

###Heatplot
d<-0.5 # try values d<1 (bimodal) or d>1 (unimodal)
B <- beta(d, d)
# density of \xi at a single point (z1, z2)
fZ <- function(z1, z2) {
  integrand <- function(c) {
    w<-(1-c^2)^(d-2)
    expo<-exp(-(z1^2+(z2-c)^2)/(2*(1-c^2)))
    w*expo
  }
  value <- integrate(integrand, lower = -1, upper = 1,
                   rel.tol = 1e-6, subdivisions = 200L)$value
  const<-1/(2*pi*2^(2*d-1)*B)
  const*value
}
# Grid for plotting in the box [-2,2]x[-2,2], modify to taste
z1_seq<-seq(-2, 2, length.out = 80)
z2_seq<-seq(-2, 2, length.out = 80)
Zmat<-outer(z1_seq, z2_seq, Vectorize(function(x, y) fZ(x, y)))
filled.contour(z1_seq, z2_seq, Zmat,
               color.palette = terrain.colors,
               xlab = "z1", ylab = "z2",
               main = bquote("Joint density of " * xi * " for " * delta == .(d)))

### 3D Plot, run previous chunk to obtain parameters
persp(z1_seq, z2_seq, Zmat,
      theta = 30, phi = 30, expand = 0.6,
      col = "blue", border="lightblue",
      xlab = "z1", ylab = "z2", zlab = "density",
      main =bquote("Joint density of " * xi * " for " * delta == .(d)))

runtime and error bound comparison for N

##### source true prob and true gradient at chunk starting in line 65
estN50<-rep(0,1000)
estN100<-rep(0,1000)
estN500<-rep(0,1000)
estN1000<-rep(0,1000)
estN5000<-rep(0,1000)
estN10000<-rep(0,1000)
estpN50<-rep(0,1000)
estpN100<-rep(0,1000)
estpN500<-rep(0,1000)
estpN1000<-rep(0,1000)
estpN5000<-rep(0,1000)
estpN10000<-rep(0,1000)
estgN50<-rep(0,1000)
estgN100<-rep(0,1000)
estgN500<-rep(0,1000)
estgN1000<-rep(0,1000)
estgN5000<-rep(0,1000)
estgN10000<-rep(0,1000)
estn50<-rep(0,1000)
estn100<-rep(0,1000)
estn500<-rep(0,1000)
estn1000<-rep(0,1000)
estn5000<-rep(0,1000)
estn10000<-rep(0,1000)
estpn50<-rep(0,1000)
estpn100<-rep(0,1000)
estpn500<-rep(0,1000)
estpn1000<-rep(0,1000)
estpn5000<-rep(0,1000)
estpn10000<-rep(0,1000)
estgn50<-rep(0,1000)
estgn100<-rep(0,1000)
estgn500<-rep(0,1000)
estgn1000<-rep(0,1000)
estgn5000<-rep(0,1000)
estgn10000<-rep(0,1000)
set.seed(123)
for(k in 1:1000){
  estradial50<-rep(0,xvalues)
  estgradient50<-rep(0,xvalues)
  estprop50<-rep(0,xvalues)
  estradial100<-rep(0,xvalues)
  estgradient100<-rep(0,xvalues)
  estprop100<-rep(0,xvalues)
  estradial500<-rep(0,xvalues)
  estgradient500<-rep(0,xvalues)
  estprop500<-rep(0,xvalues)
  estradial1000<-rep(0,xvalues)
  estgradient1000<-rep(0,xvalues)
  estprop1000<-rep(0,xvalues)
  estradial5000<-rep(0,xvalues)
  estgradient5000<-rep(0,xvalues)
  estprop5000<-rep(0,xvalues)
  estradial10000<-rep(0,xvalues)
  estgradient10000<-rep(0,xvalues)
  estprop10000<-rep(0,xvalues)
  v<-runif(10000,0,2*pi)
  c<-2*rbeta(10000,d,d)-1
  x<-seq(xlow,xup, length.out=xvalues)
  rs<-rep(0,10000)
  for(w in 1:10000){
    rs[w]<-rnorm(1,0,sqrt(1-c[w]^2))^2+rnorm(1,c[w],sqrt(1-c[w]^2))^2
  }
for(i in 1:xvalues){
  rcrit<-(-c*sin(v)+sqrt(x[i]^2+2 -c^2*cos(v)^2))/sqrt(1-c^2)
  estradial50[i]<-mean(1-exp(-rcrit[1:50]^2/2))
  estgradient50[i]<-mean(x[i]*rcrit[1:50]*exp(-rcrit[1:50]^2/2)/sqrt((x[i]^2+2 -c[1:50]^2*cos(v[1:50])^2)*(1-c[1:50]^2)))
  estprop50[i]<-sum(rs[1:50] < 2+x[i]^2)/50
  estradial100[i]<-mean(1-exp(-rcrit[1:100]^2/2))
  estgradient100[i]<-mean(x[i]*rcrit[1:100]*exp(-rcrit[1:100]^2/2)/sqrt((x[i]^2+2 -c[1:100]^2*cos(v[1:100])^2)*(1-c[1:100]^2)))
  estprop100[i]<-sum(rs[1:100] < 2+x[i]^2)/100
  estradial500[i]<-mean(1-exp(-rcrit[1:500]^2/2))
  estgradient500[i]<-mean(x[i]*rcrit[1:500]*exp(-rcrit[1:500]^2/2)/sqrt((x[i]^2+2 -c[1:500]^2*cos(v[1:500])^2)*(1-c[1:500]^2)))
  estprop500[i]<-sum(rs[1:500] < 2+x[i]^2)/500
  estradial1000[i]<-mean(1-exp(-rcrit[1:1000]^2/2))
  estgradient1000[i]<-mean(x[i]*rcrit[1:1000]*exp(-rcrit[1:1000]^2/2)/sqrt((x[i]^2+2 -c[1:1000]^2*cos(v[1:1000])^2)*(1-c[1:1000]^2)))
  estprop1000[i]<-sum(rs[1:1000] < 2+x[i]^2)/1000
  estradial5000[i]<-mean(1-exp(-rcrit[1:5000]^2/2))
  estgradient5000[i]<-mean(x[i]*rcrit[1:5000]*exp(-rcrit[1:5000]^2/2)/sqrt((x[i]^2+2 -c[1:5000]^2*cos(v[1:5000])^2)*(1-c[1:5000]^2)))
  estprop5000[i]<-sum(rs[1:5000] < 2+x[i]^2)/5000
  estradial10000[i]<-mean(1-exp(-rcrit[1:10000]^2/2))
  estgradient10000[i]<-mean(x[i]*rcrit[1:10000]*exp(-rcrit[1:10000]^2/2)/sqrt((x[i]^2+2 -c[1:10000]^2*cos(v[1:10000])^2)*(1-c[1:10000]^2)))
  estprop10000[i]<-sum(rs < 2+x[i]^2)/10000
}
  estN50[k]<-sum((trueprob-estradial50)^2)*(xup-xlow)/xvalues
  estN100[k]<-sum((trueprob-estradial100)^2)*(xup-xlow)/xvalues
  estN500[k]<-sum((trueprob-estradial500)^2)*(xup-xlow)/xvalues
  estN1000[k]<-sum((trueprob-estradial1000)^2)*(xup-xlow)/xvalues
  estN5000[k]<-sum((trueprob-estradial5000)^2)*(xup-xlow)/xvalues
  estN10000[k]<-sum((trueprob-estradial10000)^2)*(xup-xlow)/xvalues
  
  estpN50[k]<-sum((trueprob-estprop50)^2)*(xup-xlow)/xvalues
  estpN100[k]<-sum((trueprob-estprop100)^2)*(xup-xlow)/xvalues
  estpN500[k]<-sum((trueprob-estprop500)^2)*(xup-xlow)/xvalues
  estpN1000[k]<-sum((trueprob-estprop1000)^2)*(xup-xlow)/xvalues
  estpN5000[k]<-sum((trueprob-estprop5000)^2)*(xup-xlow)/xvalues
  estpN10000[k]<-sum((trueprob-estprop10000)^2)*(xup-xlow)/xvalues
  
  
  estgN50[k]<-sum((truegradient-estgradient50)^2)*(xup-xlow)/xvalues
  estgN100[k]<-sum((truegradient-estgradient100)^2)*(xup-xlow)/xvalues
  estgN500[k]<-sum((truegradient-estgradient500)^2)*(xup-xlow)/xvalues
  estgN1000[k]<-sum((truegradient-estgradient1000)^2)*(xup-xlow)/xvalues
  estgN5000[k]<-sum((truegradient-estgradient5000)^2)*(xup-xlow)/xvalues
  estgN10000[k]<-sum((truegradient-estgradient10000)^2)*(xup-xlow)/xvalues
  
  estn50[k]<-max(abs(trueprob-estradial50))
  estn100[k]<-max(abs(trueprob-estradial100))
  estn500[k]<-max(abs(trueprob-estradial500))
  estn1000[k]<-max(abs(trueprob-estradial1000))
  estn5000[k]<-max(abs(trueprob-estradial5000))
  estn10000[k]<-max(abs(trueprob-estradial10000))
  
  estpn50[k]<-max(abs(trueprob-estprop50))
  estpn100[k]<-max(abs(trueprob-estprop100))
  estpn500[k]<-max(abs(trueprob-estprop500))
  estpn1000[k]<-max(abs(trueprob-estprop1000))
  estpn5000[k]<-max(abs(trueprob-estprop5000))
  estpn10000[k]<-max(abs(trueprob-estprop10000))
  
  
  estgn50[k]<-max(abs(truegradient-estgradient50))
  estgn100[k]<-max(abs(truegradient-estgradient100))
  estgn500[k]<-max(abs(truegradient-estgradient500))
  estgn1000[k]<-max(abs(truegradient-estgradient1000))
  estgn5000[k]<-max(abs(truegradient-estgradient5000))
  estgn10000[k]<-max(abs(truegradient-estgradient10000))
}

par(mfrow=c(2,2))
boxplot(estN50,estN100,estN500,estN1000,estN5000,estN10000, names=c("50","100","500","1000","5000","10000"), main=bquote("Integrated Squared Error of "*Phi[N]), ylim=c(0,0.006), col = "red1")
boxplot(estpN50,estpN100,estpN500,estpN1000,estpN5000,estpN10000, names=c("50","100","500","1000","5000","10000"), main=bquote("Integrated Squared Error of "*tilde(Phi[N])), ylim=c(0,0.006), col = "red1")

boxplot(log(estN50),log(estN100),log(estN500),log(estN1000),log(estN5000),log(estN10000), names=c("50","100","500","1000","5000","10000"), main=bquote("(Log)-Integrated Squared Error of "*Phi[N]),col="palevioletred", ylim=c(-23,-3))
boxplot(log(estpN50),log(estpN100),log(estpN500),log(estpN1000),log(estpN5000),log(estpN10000), names=c("50","100","500","1000","5000","10000"), main=bquote("(Log)-Integrated Squared Error of "*tilde(Phi)[N]),col="palevioletred",ylim=c(-23,-3))

par(mfrow=c(2,1))
boxplot(estgN50,estgN100,estgN500,estgN1000,estgN5000,estgN10000, names=c("50","100","500","1000","5000","10000"), main=bquote("Integrated Squared Error of "*nabla*Phi[N]), col = "plum1")
boxplot(log(estgN50),log(estgN100),log(estgN500),log(estgN1000),log(estgN5000),log(estgN10000), names=c("50","100","500","1000","5000","10000"), main=bquote("(Log-) Integrated Squared Error of "*nabla*Phi[N]), col = "plum1")

par(mfrow=c(2,2))
boxplot(estn50,estn100,estn500,estn1000,estn5000,estn10000, names=c("50","100","500","1000","5000","10000"), main=bquote("Maximum Error of "*Phi[N]), ylim=c(0,0.1), col = "blue")
boxplot(estpn50,estpn100,estpn500,estpn1000,estpn5000,estpn10000, names=c("50","100","500","1000","5000","10000"), main=bquote("Maximum Error of "*tilde(Phi[N])), ylim=c(0,0.1), col = "blue")

boxplot(log(estn50),log(estn100),log(estn500),log(estn1000),log(estn5000),log(estn10000), names=c("50","100","500","1000","5000","10000"), main=bquote("(Log)-Maximum Error of "*Phi[N]),col="lightblue3", ylim=c(-23,-3))
boxplot(log(estpn50),log(estpn100),log(estpn500),log(estpn1000),log(estpn5000),log(estpn10000), names=c("50","100","500","1000","5000","10000"), main=bquote("(Log)-Maximum Error of "*tilde(Phi)[N]),col="lightblue3",ylim=c(-23,-3))

par(mfrow=c(2,1))
boxplot(estgn50,estgn100,estgn500,estgn1000,estgn5000,estgn10000, names=c("50","100","500","1000","5000","10000"), main=bquote("Integrated Squared Error of "*nabla*Phi[N]), col = "skyblue1")
boxplot(log(estgn50),log(estgn100),log(estgn500),log(estgn1000),log(estgn5000),log(estgn10000), names=c("50","100","500","1000","5000","10000"), main=bquote("(Log-) Integrated Squared Error of "*nabla*Phi[N]), col = "skyblue1")