inventory_simul<-function(lostsalescosts=10,holdingcosts=0.1,upto=65){
#returns the total inventory costs

ndays<-10;
demand<-rgamma(ndays,5);
startiv <- 25;
dayorderarrives <- 5;
lostsales<-rep(0,10);
ei<-startiv-demand[1];
for(i in 2:10){
  ei[i]<- ei[i-1]-demand[i];
  if(i==dayorderarrives) ei[i] = ei[i]+upto-startiv;
  if(ei[i]<0){
	 lostsales[i]= -ei[i];
	  ei[i] =0;
  }
}
costs <- sum(ei)*holdingcosts+sum(lostsales)*lostsalescosts;
#print(dayorderarrives);
#print(data.frame(demand,ei,lostsales));

costs
}
#inventory_simul(10,0.1,65)

vect_simul_naive<-function(n,...){
#returns mean and se of the vector of total inventory costs
res<-1;
for(i in  1:n){res[i]<- inventory_simul(...)}
res1<-mean(res);
res1[2]<-sd(res)/sqrt(n);
res1
}


qsimul <- function(demand,lostsalescosts=10,holdingcosts=0.1,upto=65){
#evaluates the functionq(..) as function of the random input

ndays<-10;
startiv <- 25;
lostsales<-rep(0,10);
dayorderarrives <- 5;
ei<-startiv-demand[1];
for(i in 2:10){
  ei[i]<- ei[i-1]-demand[i];
  if(i==dayorderarrives) ei[i] = ei[i]+upto-startiv;
  if(ei[i]<0){
	 lostsales[i]= -ei[i];
	  ei[i] =0;
  }
}

costs <- sum(ei)*holdingcosts+sum(lostsales)*lostsalescosts;
#print(dayorderarrives);
#print(data.frame(demand,ei,lostsales));

costs


}


inventory_simul_IS<-function(betais,...){
#returns the total inventory costs
# IS with gamma(5) distribution with beta < 1

demand<-rgamma(ndays,5,betais);

qx <- qsimul(demand,...);

fx<-dgamma(demand,5);
gx<-dgamma(demand,5,betais);
w<-prod(fx/gx);

w*qx
}

vect_simul_IS<-function(n,betais,...){
#returns mean and se of the vector of total inventory costs
res<-1;
for(i in  1:n){res[i]<- inventory_simul_IS(betais,...)}
res1<-mean(res);
res1[2]<-sd(res)/sqrt(n);
res1
}



> vect_simul_IS(10000,1.01)
[1] 25.6412820  0.1390444
> vect_simul_IS(10000,.98)
[1] 25.9133288  0.1150556
> vect_simul_IS(10000,.98)
[1] 25.6878865  0.1106122
> vect_simul_IS(10000,1.)
[1] 25.7531443  0.1263299
> vect_simul_IS(10000,1.)
[1] 25.8433231  0.1303233
> vect_simul_IS(10000,.96)
[1] 25.8435170  0.1177698
> vect_simul_IS(10000,.96)
[1] 25.8930192  0.1187817
> vect_simul_IS(10000,.96)
[1] 25.8663792  0.1177927
> vect_simul_IS(10000,.97)
[1] 25.8191411  0.1130315
> vect_simul_IS(10000,.97)
[1] 26.0434391  0.1148278
> vect_simul_IS(10000,.97)
[1] 25.8950951  0.1122108
> (0.13/0.113)^2
[1] 1.323518