set.seed(12345)

n<-10000;

z1<-rnorm(n,0,1)
z2<-rnorm(n,0,1)
z3<-rnorm(n,0,1)
z4<-rnorm(n,0,1)
z5<-rnorm(n,0,1)
z6<-rnorm(n,0,1)
z7<-rnorm(n,0,1)
z8<-rnorm(n,0,1)
z9<-rnorm(n,0,1)
z10<-rnorm(n,0,1)

chi15<-z1^2+z2^2+z3^2+z4^2+z5^2
chi25<-z6^2+z7^2+z8^2+z9^2+z10^2

f55<-(chi15/5)/(chi25/5)

hist(f55, breaks=100, xlim=c(0, 20))

par(new=T)
plot(density(f55), axes=F, main="", xlim=c(0, 20), lwd=2, col="blue")

df1_list<-c(5,9,15,38)
df2_list<-c(10,10,20,40)
par(mfrow=c(2,2))

for (i in 1:length(df1_list)) {
  hist(rf(n, df1=df1_list[i], df2=df2_list[i], ncp=0), breaks=100, xlab="f", main=paste("df1=", df1_list[i], "df2=", df2_list[i]))
}

par(mfrow=c(1,1))
curve(df(x, 3, 15, ncp=0), col="blue", xlim=c(0,5), ylim=c(0,1), xlab="f", ylab="f(f)")
curve(df(x, 5, 15, ncp=0), add=T, col="red", xlim=c(0,5), ylim=c(0,1), xlab="f", ylab="f(f)")
curve(df(x, 10, 15, ncp=0), add=T, col="green", xlim=c(0,5), ylim=c(0,1), xlab="f", ylab="f(f)")
curve(df(x, 15, 15, ncp=0), add=T, col="black", xlim=c(0,5), ylim=c(0,1), xlab="f", ylab="f(f)")

f11<-rep(NA,10)
f12<-rep(NA,10)
f13<-rep(NA,10)
f14<-rep(NA,10)
f15<-rep(NA,10)
f16<-rep(NA,10)
f17<-rep(NA,10)
f18<-rep(NA,10)
f19<-rep(NA,10)
f110<-rep(NA,10)

for(i in 1:10) {f11[i]<-qf(0.95, 1, i)}
for(i in 1:10) {f12[i]<-qf(0.95, 2, i)}
for(i in 1:10) {f13[i]<-qf(0.95, 3, i)}
for(i in 1:10) {f14[i]<-qf(0.95, 4, i)}
for(i in 1:10) {f15[i]<-qf(0.95, 5, i)}
for(i in 1:10) {f16[i]<-qf(0.95, 6, i)}
for(i in 1:10) {f17[i]<-qf(0.95, 7, i)}
for(i in 1:10) {f18[i]<-qf(0.95, 8, i)}
for(i in 1:10) {f19[i]<-qf(0.95, 9, i)}
for(i in 1:10) {f110[i]<-qf(0.95, 10, i)}

round((poi<-cbind(f11,f12,f13,f14,f15,f16,f17,f18,f19,f110)), digits=2)