#The following function takes a number of samples (n, default: 10) and a number of degrees of freedom (df, default:2), and generates n samples from a distribution that is the natural logarithm of the sum of df standard normal distributions.

rJimenez<-function(n=10, df=2){
  y<-replicate(n,{
    x=rnorm(df,0,1)
    y=sum(log(abs(x)))
  })
  return(y)
}

#We generate 100000 samples from a Jimenez distribution with 10 degrees of freedom:

Jimenez<-rJimenez(100000,10)

#We graph a series of plots from our sample of Jimenez distributions with 10 degrees of freedom

par(mfrow=c(2,2), mar=c(3,4,1,1))#the matrix of plots is set up
plot(Jimenez)#a scatter plot
hist(Jimenez)#a histogram
plot(density(Jimenez))#a density plot
plot(ecdf(Jimenez))#a cumulative density plot

#The following is a function that samples n (default: 10) F distributions with df1 (default: 1) numerator degrees of freedom and df2 (default: 2) denominator degrees of freedom.

rfJimenez<-function(n=10, df1=1, df2=2){
  w<-replicate(n,{
    x=rnorm(df1,0,1)
    chi1=sum(x^2)
    z=rnorm(df2,0,1)
    chi2=sum(z^2)
    y=(chi1/df1)/(chi2/df2)
  })
  return(w)
}

#We now use the rfJimenez function to generate samples of 10, 100, 1000, and 10000 from F distributions with 10 numerator degrees of freedom and 20 denominator degrees of freedom

x10<-rfJimenez(10,10,20)

x100<-rfJimenez(100,10,20)

x1000<-rfJimenez(1000,10,20)

x100000<-rfJimenez(100000,10,20)

#The following function takes a sample from an F distribution and empirically test the hypothesis that the sample's mean corresponds with the theoretical mean

Mtest<-function(x, n=10, k=10000, df1=10, df2=20){
  xt<-mean(x)-(df2/(df2-2))
  y=replicate(k,{y=mean(rf(n,df1,df2))-(df2/(df2-2))})
  P=length(y[y>xt])/k
  return(P)
}

#We now use the Mtest function to test the samples generated with the rfJimenez function

Mtest(x10)

Mtest(x100,n=100)

Mtest(x1000,n=1000)

Mtest(x100000,n=100000)

#The following function takes a sample from an F distribution and empirically test the hypothesis that the sample's variance corresponds with the theoretical variance

Vtest<-function(x, n=10, k=10000, df1=10, df2=20){
  V=((2*(df2^2)*(df1+df2-2))/(df1*((df2-2)^2)*(df2-4)))
  xt<-var(x)-V
  y=replicate(k,{y=var(rf(n,df1,df2))-V})
  P=length(y[y>xt])/k
  return(P)
}

#We now use the Vtest function to test the samples generated with the rfJimenez function

Vtest(x10)

Vtest(x100,n=100)

Vtest(x1000,n=1000)

Vtest(x100000,n=100000)