Monte Carlo Optimization, Simulation and Sensitivity of Queueing Networks,  Reuven Y.Rubinstein, Wiley

```{r}

#p.72

a=20

g=function(y){

return(exp(-y^2)-a)

}


x=0.01
ite=10^4
h=0.001

for(l in 1:ite){

dg=(g(x+h)-g(x))/h

x=x-eta*(g(x)-a)/dg

print(g(x))

}



```


```{r}


#p.73

#Robin and Monro(1951)

a=20

g=function(y){

return(exp(-y^2))

}

pthi=function(z){

return(g(z)+a+rnorm(1,0,1))

}

#a_i=1/i,i=1,2,.....

x=0.01
ite=10^4

sita=198
ep=0.1

vectors=c()

for(l in 1:ite){

x_pre=x

x=x-(pthi(x))*(1/l)

vectors=c(vectors,abs(x_pre-x))

R=mean(pthi(rep(x,10000)))



}

R=mean(pthi(rep(x,10000)))

#Theorem 1

if(abs(abs(x)-sita)<1/ep & abs(abs(x)-sita)>ep){

print((abs(x)-sita)*(R-a))
  
}  

#d
mean(pthi(rep(x,10000))^2)/(1+x^2)



```

```{r}

#p.122
#Normal

pthi=function(z){
  
return(exp(-(z-2)^2))  
  
}

mu=1
sig=1
N=10000

ite=1000
eta=10^(-2)

for(l in 1:ite){
  
Y=rnorm(N,mu,sig^2)  
  
dgmu=(Y-mu)/(sig^2)  
dgsig=-log(sig)-(Y-mu)^2/(sig^2)

mu=mu+eta*mean(pthi(Y)*dgmu)
sig=sig+eta*mean(pthi(Y)*dgsig)
  
g0=mean(pthi(Y)*dnorm(Y,mu,sig^2))

print(g0)

}


```

```{r}

#p.122
#exponential

pthi=function(z){
  
return(exp(-(z-2)^2))  
  
}

lam=0.01
N=100000

ite=10^3
eta=10^(-2)

for(l in 1:ite){
  
Y=rexp(N,lam)  
  
dlam=(1/lam-Y)
ddlam=-1/(lam^2)+dlam^2
lam=lam-eta*mean(pthi(Y)*dlam)/mean(pthi(Y)*ddlam)
g0=mean(pthi(Y)*dexp(Y,lam))

print(g0)

}


```

```{r}

#p.122
#Weibull

pthi=function(z){
  
return(exp(-(z-2)^2))  
  
}

alpha=1
beta=2
N=100000

ite=10^3
eta=10^(-1)

for(l in 1:ite){
  
Y=rweibull(N,alpha,beta)  
  
dalpha=1/alpha-log(beta)+log(Y)-((Y/beta)^alpha)*log(Y/beta)
dbeta=-alpha/beta+(alpha/beta)*(Y/beta)^alpha

alpha=alpha+eta*mean(pthi(Y)*dalpha)
beta=beta+eta*mean(pthi(Y)*dbeta)

g0=mean(pthi(Y)*dweibull(Y,alpha,beta))

print(g0)

}


```

```{r}

#p.206

#RS-1(the random search double-trial algorithm)

n=100
p=10
X=matrix(runif(n*p,-2,2),nrow=n,ncol=p)
b=c(2,6,-3,1,4,rep(0,5))
y=X%*%b+rnorm(n,0,1)

g=function(B){
return(sum((y-X%*%B)^2))  
}

alpha=0.01
beta=0.001
coef=rep(0.1,p)
ite=100

for(l in 1:ite){
  
for(i in 1:p){
x1=coef;x2=coef
sig=runif(1,0,1)
x1[i]=x1[i]+beta*sig
x2[i]=x2[i]-beta*sig
coef[i]=coef[i]-alpha*(g(x1)-g(x2))*sig/(2*beta)
}

print(g(coef))
 
}


```


```{r}

#p.206

#RS-2(acceptance-rejection random search algorithm)

n=100
p=10
X=matrix(runif(n*p,-2,2),nrow=n,ncol=p)
b=c(2,6,-3,1,4,rep(0,5))
y=X%*%b+rnorm(n,0,1)

g=function(B){
return(sum((y-X%*%B)^2))  
}

alpha=0.01
beta=0.001
coef=rep(0.1,p)
ite=100

for(l in 1:ite){
  
for(i in 1:p){
x1=coef;x2=coef
sig=runif(1,0,1)
x1[i]=x1[i]+beta*sig
if(g(x1)-g(x2)<0){
coef[i]=coef[i]-alpha*(g(x1)-g(x2))*sig/beta
}
}
print(g(coef))
}


```


```{r}

#p.206

#RS-3(optimal trial random search algorithm)

n=10
p=3
X=matrix(runif(n*p,-2,2),nrow=n,ncol=p)
b=c(2,6,-3,1,4,rep(0,5))
b=b[1:p]
y=X%*%b+rnorm(n,0,1)

g=function(B){
return(sum((y-X%*%B)^2))  
}

alpha=0.1
coef=rep(0.1,p)
N=100
ite=2*10^2

for(l in 1:ite){
for(i in 1:p){
beta=sample(seq(0,1,0.1),N,replace=T)
sig=runif(N,0,1)  
yvec=c()
for(k in 1:N){
x1=coef;x2=coef
x1[i]=x1[i]+beta[k]
yvec=c(yvec,g(x1)-g(x2))
}
no=(c(1:N)[yvec==min(yvec)])[1]
coef[i]=coef[i]-alpha*beta[no]*(yvec[no])*sig[no]
}
print(g(coef))
 
}


```


```{r}

#p.206

#RS-4(statistical gradient random search algorithm)

n=10
p=3
X=matrix(runif(n*p,-2,2),nrow=n,ncol=p)
b=c(2,6,-3,1,4,rep(0,5))
b=b[1:p]
y=X%*%b+rnorm(n,0,1)

g=function(B){
return(sum((y-X%*%B)^2))  
}

alpha=0.01
coef=rep(0.1,p)
N=1000
ite=3*10^3

for(l in 1:ite){
for(i in 1:p){
beta=sample(seq(0,1,0.1),N,replace=T)
sig=runif(N,0,1)  
yvec=c()
for(k in 1:N){
x1=coef;x2=coef
x1[i]=x1[i]+beta[k]
yvec=c(yvec,g(x1)-g(x2))
}
beta=sample(1:100,1)
coef[i]=coef[i]-alpha*(1/beta)*mean(yvec*sig)
}
print(g(coef))
 
}


```