# R/rMHN_negative_gamma.R In subhadippal2019/MHN: Generates Random Variables from the Modified Half Normal Distribution

#### Defines functions rMHN_negative_b_altrMHN_negative_b

```  #samplerG for negative x values

#' @export
#'
rMHN_negative_b<-function(n=1, alpha,beta,gamma){

a=beta; b=gamma
# b=abs(b)
if(1-(a>0)*(b<=0)*(alpha>0)){  print("The parameters beta , alpha have to be positive and gamma has to be nonpositive to use this functioin."); return(NULL) }
cut_off_variable=1
# cut_off_variable: decides which functional approximations to be used to crate the upper bound for the rejection sampler.
b=abs(b)

if(alpha<=cut_off_variable){ # Applly Strategy1
shape_par=(a+b)/(2*a+b)*alpha; rate_par=a+b; Exponent_power=(a+b)/(2*a+b);
Exp_g<-(  rgamma(n = n,shape=shape_par,rate=rate_par) ) ##may use rgamss
g=Exp_g^(Exponent_power)

log_U=log(runif(n = n, min = 0, max = 1))
log_Rejection_Rate = (  -a*g^2 -  b*g +  (a+b)*Exp_g     )
sample_x <- g[log_U<=log_Rejection_Rate]

size=n-length(sample_x)
i=1
while(size!=0){
Exp_g<-(  rgamma(n = size,shape=shape_par,rate=rate_par) ) ##may use rgamss
g=Exp_g^(Exponent_power)

log_U=log(runif(n = size, min = 0, max = 1))
log_Rejection_Rate = (  -a*g^2 -  b*g +  (a+b)*Exp_g     )
sample_x=c(sample_x,g[log_U<=log_Rejection_Rate])
size=n-length(sample_x)
i=i+1
#print(i)
}

x=sample_x
}

if(alpha>cut_off_variable){ # Applly Strategy1
#x^(alpha-1)exp(-ax^2-bx). a>0,b>0, alpha>1
x_mode=(sqrt(b^2+8*(alpha-1)*a)-b)/(4*a)
scale=x_mode
#Scaled_rv= original_rv/scale; original_rv=Scaled_rv*scale
new_a=a* (scale)^2;new_b=b*(scale);
shape_par=(new_a+new_b)/(2*new_a+new_b)*alpha; rate_par=new_a+new_b; Exponent_power=(new_a+new_b)/(2*new_a+new_b);

Exp_g<-(  rgamma(n = n,shape=shape_par,rate=rate_par) ) ##may use rgamss
g=Exp_g^(Exponent_power)

log_U=log(runif(n = n, min = 0, max = 1))
log_Rejection_Rate =  (  -new_a*g^2 -  new_b*g +  (new_a+new_b)*Exp_g)
sample_x <- g[log_U<=log_Rejection_Rate]
size=n-length(sample_x)
i=1
while(size!=0){
Exp_g<-(  rgamma(n = size,shape=shape_par,rate=rate_par) ) ##may use rgamss
g=Exp_g^(Exponent_power)

log_U=log(runif(n = size, min = 0, max = 1))
log_Rejection_Rate =  (  -new_a*g^2 -  new_b*g +  (new_a+new_b)*Exp_g)
#sample_x <- x[log_U<=log_Rejection_Rate]
sample_x=c(sample_x,g[log_U<=log_Rejection_Rate])
size=n-length(sample_x)
i=i+1
#print(i)
}

x=sample_x*scale
}
return(x)
}

rMHN_negative_b_alt<-function(alpha,beta,gamma){

a=beta; b=gamma
# b=abs(b)
if(1-(a>0)*(b<=0)*(alpha>0)){  print("The parameters beta , alpha have to be positive and gamma has to be nonpositive to use this functioin."); return(NULL) }
cut_off_variable=1
# cut_off_variable: decides which functional approximations to be used to crate the upper bound for the rejection sampler.
b=abs(b)
ACCEPT=FALSE

if(alpha<=cut_off_variable){ # Applly Strategy1
#print("here")
counter=0;
while(!ACCEPT){
counter=counter+1
shape_par=(a+b)/(2*a+b)*alpha; rate_par=a+b; Exponent_power=(a+b)/(2*a+b);

Exp_g<-(  rgamma(1,shape=shape_par,rate=rate_par) ) ##may use rgamss
g=Exp_g^(Exponent_power)
if(     log(runif(1)) <=  (  -a*g^2 -  b*g +  (a+b)*Exp_g     )  ){
ACCEPT=T
#print(paste("STEP1=",counter))
return(g)
}
if(counter>100000){print("NO acceptence after 100000 iteration Procedure1");ACCEPT=T;}
}
}

if(alpha>cut_off_variable){ # Applly Strategy1
#x^(alpha-1)exp(-ax^2-bx). a>0,b>0, alpha>1
x_mode=(sqrt(b^2+8*(alpha-1)*a)-b)/(4*a)
scale=x_mode
#Scaled_rv= original_rv/scale; original_rv=Scaled_rv*scale
new_a=a* (scale)^2;new_b=b*(scale);
ACCEPT=F
counter=0;
while(!ACCEPT){
counter=counter+1
shape_par=(new_a+new_b)/(2*new_a+new_b)*alpha; rate_par=new_a+new_b; Exponent_power=(new_a+new_b)/(2*new_a+new_b);

Exp_g<-(  rgamma(1,shape=shape_par,rate=rate_par) ) ##may use rgamss
g=Exp_g^(Exponent_power)
if(     log(runif(1)) <=  (  -new_a*g^2 -  new_b*g +  (new_a+new_b)*Exp_g)  ){
ACCEPT=T
# print(paste("STEP2=",counter))
x_sampled=g*scale
return(x_sampled)
}
if(counter>100000){print("NO acceptence after 100000 iteration Procedure1");ACCEPT=T;}
}
}
return(NULL)
}
```
subhadippal2019/MHN documentation built on Jan. 10, 2022, 12:11 p.m.