c ***************************************************************** c Summarize the simulation for method of moments c ***************************************************************** subroutine momsum(ansmom,dummysim,Nbigsim,smsumer) implicit real*8(a-h,o-z) real*8 ansmom(500,5),dummysim(500),smsumer(2,7) print *, " " print *, "Results of the simulation, Moments" print *, "beta0 = ",(ansmom(i,3),i=1,Nbigsim) b0mean = smean(ansmom,500,5,3,1,Nbigsim) b0sd = stmean(ansmom,500,5,3,1,Nbigsim) b0med = smedian(ansmom,dummysim,500,5,3,1,Nbigsim) b0mad = smad(ansmom,dummysim,500,5,3,1,Nbigsim) print *, "The mean of the beta0 = ",b0mean print *, "The s.d. of the beta0 = ",b0sd print *, "The median of the beta0 = ",b0med print *, "The mad of the beta0 = ",b0mad print *, " " print *, "Results of the simulation, Moments" print *, "beta1 = ",(ansmom(i,4),i=1,Nbigsim) b1mean = smean(ansmom,500,5,4,1,Nbigsim) b1sd = stmean(ansmom,500,5,4,1,Nbigsim) b1med = smedian(ansmom,dummysim,500,5,4,1,Nbigsim) b1mad = smad(ansmom,dummysim,500,5,4,1,Nbigsim) b1mae = smae(ansmom,dummysim,500,5,4,1,Nbigsim,0.83d0) b1mse = ((b1mean - 0.83) ** 2) + (b1sd **2) b1ase = ((b1med - 0.83) ** 2) + (b1mad **2) print *, "The mean of the beta1 = ",b1mean print *, "The s.d. of the beta1 = ",b1sd print *, "The mse of the beta1 = ",b1mse print *, "The mae of the beta1 = ",b1mae print *, "The median of the beta1 = ",b1med print *, "The mad of the beta1 = ",b1mad print *, "The madmse of the beta1 = ",b1ase print *, " " print *, "Results of the simulation, Moments" print *, "sigma_e^2 = ",(ansmom(i,1),i=1,Nbigsim) semean = smean(ansmom,500,5,1,1,Nbigsim) sesd = stmean(ansmom,500,5,1,1,Nbigsim) semed = smedian(ansmom,dummysim,500,5,1,1,Nbigsim) semad = smad(ansmom,dummysim,500,5,1,1,Nbigsim) print *, "The mean of the sigma_e^2 = ",semean print *, "The s.d. of the sigma_e^2 = ",sesd print *, "The median of the sigma_e^2 = ",semed print *, "The mad of the sigma_e^2 = ",semad print *, " " print *, "Results of the simulation, Moments" print *, "sigma_u^2 = ",(ansmom(i,2),i=1,Nbigsim) sumean = smean(ansmom,500,5,2,1,Nbigsim) susd = stmean(ansmom,500,5,2,1,Nbigsim) sumed = smedian(ansmom,dummysim,500,5,2,1,Nbigsim) sumad = smad(ansmom,dummysim,500,5,2,1,Nbigsim) print *, "The mean of the sigma_u^2 = ",sumean print *, "The s.d. of the sigma_u^2 = ",susd print *, "The median of the sigma_u^2 = ",sumed print *, "The mad of the sigma_u^2 = ",sumad print *, " " print *, "Results of the simulation, Moments" print *, "rhoqt = ",(ansmom(i,5),i=1,Nbigsim) rhomean = smean(ansmom,500,5,5,1,Nbigsim) rhosd = stmean(ansmom,500,5,5,1,Nbigsim) rhomed = smedian(ansmom,dummysim,500,5,5,1,Nbigsim) rhomad = smad(ansmom,dummysim,500,5,5,1,Nbigsim) rhomae = smae(ansmom,dummysim,500,5,5,1,Nbigsim,0.54d0) rhomse = ((rhomean - 0.54) ** 2) + (rhosd **2) rhoase = ((rhomed - 0.54) ** 2) + (rhomad **2) print *, "The mean of the rho = ",rhomean print *, "The s.d. of the rho = ",rhosd print *, "The mse of the rho = ",rhomse print *, "The mae of the rho = ",rhomae print *, "The median of the rho = ",rhomed print *, "The mad of the rho = ",rhomad print *, "The madmse of the rho = ",rhoase do 455 i=1,500 cc = ansmom(i,5) if (dabs(cc).ge.0.90d0) cc = 0.90 cc = dlog((1.0d0 + cc) / (1.0d0 - cc)) 455 ansmom(i,5) = cc cc = dlog((1.0d0 + 0.54d0) / (1.0d0 - 0.54d0)) print *, " " print *, "Results of the simulation, Moments" print *, "Fisher rhoqt = ",(ansmom(i,5),i=1,Nbigsim) rrhomean = smean(ansmom,500,5,5,1,Nbigsim) rrhosd = stmean(ansmom,500,5,5,1,Nbigsim) rrhomed = smedian(ansmom,dummysim,500,5,5,1,Nbigsim) rrhomad = smad(ansmom,dummysim,500,5,5,1,Nbigsim) rrhomae = smae(ansmom,dummysim,500,5,5,1,Nbigsim,cc) rrhomse = ((rrhomean - cc) ** 2) + (rrhosd **2) rrhoase = ((rrhomed - cc) ** 2) + (rrhomad **2) print *, "The mean of the Fisher rho = ",rrhomean print *, "The s.d. of the Fisher rho = ",rrhosd print *, "The mse of the Fisher rho = ",rrhomse print *, "The mae of the Fisher rho = ",rrhomae print *, "The median of the Fisher rho = ",rrhomed print *, "The mad of the Fisher rho = ",rrhomad print *, "The madmse of the Fisher rho = ",rrhoase smsumer(1,1) = b1mean smsumer(1,2) = b1sd smsumer(1,3) = b1med smsumer(1,4) = b1mad smsumer(1,5) = b1mae smsumer(1,6) = b1mse smsumer(1,7) = b1ase smsumer(2,1) = rhomean smsumer(2,2) = rhosd smsumer(2,3) = rhomed smsumer(2,4) = rhomad smsumer(2,5) = rhomae smsumer(2,6) = rhomse smsumer(2,7) = rhoase return end c ***************************************************************** c Summarize the simulation for Bayes c ***************************************************************** subroutine simsum(anssim,dummysim,Nbigsim,smsumer) implicit real*8(a-h,o-z) real*8 anssim(500,6),dummysim(500),smsumer(2,7) print *, " " print *, "Results of the simulation, Bayes" print *, "beta0 = ",(anssim(i,3),i=1,Nbigsim) b0mean = smean(anssim,500,6,3,1,Nbigsim) b0sd = stmean(anssim,500,6,3,1,Nbigsim) b0med = smedian(anssim,dummysim,500,6,3,1,Nbigsim) b0mad = smad(anssim,dummysim,500,6,3,1,Nbigsim) print *, "The mean of the beta0 = ",b0mean print *, "The s.d. of the beta0 = ",b0sd print *, "The median of the beta0 = ",b0med print *, "The mad of the beta0 = ",b0mad print *, " " print *, "Results of the simulation, Bayes" print *, "beta1 = ",(anssim(i,4),i=1,Nbigsim) b1mean = smean(anssim,500,6,4,1,Nbigsim) b1sd = stmean(anssim,500,6,4,1,Nbigsim) b1med = smedian(anssim,dummysim,500,6,4,1,Nbigsim) b1mad = smad(anssim,dummysim,500,6,4,1,Nbigsim) b1mae = smae(anssim,dummysim,500,6,4,1,Nbigsim,0.83d0) b1mse = ((b1mean - 0.83) ** 2) + (b1sd **2) b1ase = ((b1med - 0.83) ** 2) + (b1mad **2) print *, "The mean of the beta1 = ",b1mean print *, "The s.d. of the beta1 = ",b1sd print *, "The mse of the beta1 = ",b1mse print *, "The mae of the beta1 = ",b1mae print *, "The median of the beta1 = ",b1med print *, "The mad of the beta1 = ",b1mad print *, "The madmse of the beta1 = ",b1ase print *, " " print *, "Results of the simulation, Bayes" print *, "sigma_e^2 = ",(anssim(i,1),i=1,Nbigsim) semean = smean(anssim,500,6,1,1,Nbigsim) sesd = stmean(anssim,500,6,1,1,Nbigsim) semed = smedian(anssim,dummysim,500,6,1,1,Nbigsim) semad = smad(anssim,dummysim,500,6,1,1,Nbigsim) print *, "The mean of the sigma_e^2 = ",semean print *, "The s.d. of the sigma_e^2 = ",sesd print *, "The median of the sigma_e^2 = ",semed print *, "The mad of the sigma_e^2 = ",semad print *, " " print *, "Results of the simulation, Bayes" print *, "sigma_u^2 = ",(anssim(i,2),i=1,Nbigsim) sumean = smean(anssim,500,6,2,1,Nbigsim) susd = stmean(anssim,500,6,2,1,Nbigsim) sumed = smedian(anssim,dummysim,500,6,2,1,Nbigsim) sumad = smad(anssim,dummysim,500,6,2,1,Nbigsim) print *, "The mean of the sigma_u^2 = ",sumean print *, "The s.d. of the sigma_u^2 = ",susd print *, "The median of the sigma_u^2 = ",sumed print *, "The mad of the sigma_u^2 = ",sumad print *, " " print *, "Results of the simulation, Bayes" print *, "rhoqt = ",(anssim(i,5),i=1,Nbigsim) rhomean = smean(anssim,500,6,5,1,Nbigsim) rhosd = stmean(anssim,500,6,5,1,Nbigsim) rhomed = smedian(anssim,dummysim,500,6,5,1,Nbigsim) rhomad = smad(anssim,dummysim,500,6,5,1,Nbigsim) print *, "The mean of the rho = ",rhomean print *, "The s.d. of the rho = ",rhosd print *, "The median of the rho = ",rhomed print *, "The mad of the rho = ",rhomad print *, " " print *, "Results of the simulation, Bayes" print *, "rhoqt = ",(anssim(i,6),i=1,Nbigsim) rhomean = smean(anssim,500,6,6,1,Nbigsim) rhosd = stmean(anssim,500,6,6,1,Nbigsim) rhomed = smedian(anssim,dummysim,500,6,6,1,Nbigsim) rhomad = smad(anssim,dummysim,500,6,6,1,Nbigsim) rhomae = smae(anssim,dummysim,500,6,6,1,Nbigsim,0.54d0) rhomse = ((rhomean - 0.54) ** 2) + (rhosd **2) rhoase = ((rhomed - 0.54) ** 2) + (rhomad **2) print *, "The mean of the rhonew = ",rhomean print *, "The s.d. of the rhonew = ",rhosd print *, "The mse of the rhonew = ",rhomse print *, "The mae of the rhonew = ",rhomae print *, "The median of the rhonew = ",rhomed print *, "The mad of the rhonew = ",rhomad print *, "The madmse of the rhonew = ",rhoase smsumer(1,1) = b1mean smsumer(1,2) = b1sd smsumer(1,3) = b1med smsumer(1,4) = b1mad smsumer(1,5) = b1mae smsumer(1,6) = b1mse smsumer(1,7) = b1ase smsumer(2,1) = rhomean smsumer(2,2) = rhosd smsumer(2,3) = rhomed smsumer(2,4) = rhomad smsumer(2,5) = rhomae smsumer(2,6) = rhomse smsumer(2,7) = rhoase do 455 i=1,500 cc = anssim(i,6) if (dabs(cc).ge.0.90d0) cc = 0.90 cc = dlog((1.0d0 + cc) / (1.0d0 - cc)) 455 anssim(i,6) = cc cc = dlog((1.0d0 + 0.54d0) / (1.0d0 - 0.54d0)) print *, " " print *, "Results of the simulation, Bayes" print *, "Fisher rhoqt = ",(anssim(i,6),i=1,Nbigsim) rrhomean = smean(anssim,500,6,6,1,Nbigsim) rrhosd = stmean(anssim,500,6,6,1,Nbigsim) rrhomed = smedian(anssim,dummysim,500,6,6,1,Nbigsim) rrhomad = smad(anssim,dummysim,500,6,6,1,Nbigsim) rrhomae = smae(anssim,dummysim,500,6,6,1,Nbigsim,cc) rrhomse = ((rrhomean - cc) ** 2) + (rrhosd **2) rrhoase = ((rrhomed - cc) ** 2) + (rrhomad **2) print *, "The mean of the Fisher rho = ",rrhomean print *, "The s.d. of the Fisher rho = ",rrhosd print *, "The mse of the Fisher rho = ",rrhomse print *, "The mae of the Fisher rho = ",rrhomae print *, "The median of the Fisher rho = ",rrhomed print *, "The mad of the Fisher rho = ",rrhomad print *, "The madmse of the Fisher rho = ",rrhoase return end c ***************************************************************** c Summarize the simulation for Bayes, number of mixtures c ***************************************************************** subroutine simpop(npopsim,Nbigsim) implicit real*8(a-h,o-z) integer npopsim(500,1) real*8 ppop(5) do 100 j=1,5 100 ppop(j) = 0. do 200 i=1,Nbigsim do 150 j=1,5 if (npopsim(i,1).eq.j) ppop(j) = ppop(j) + (1.0d0 / Nbigsim) 150 continue 200 continue do 300 j=1,5 print *, "Kbigest and % it was selected = ", 6 j,ppop(j) 300 continue return end c ***************************************************************** c Summarize the simulation for True MLE c ***************************************************************** subroutine mletrsum(ansmletr,dummysim,Nbigsim,smsumer) implicit real*8(a-h,o-z) real*8 ansmletr(500,5),dummysim(500),smsumer(2,7) print *, " " print *, "Results of the simulation, True MLE" print *, "beta0 = ",(ansmletr(i,3),i=1,Nbigsim) b0mean = smean(ansmletr,500,5,3,1,Nbigsim) b0sd = stmean(ansmletr,500,5,3,1,Nbigsim) b0med = smedian(ansmletr,dummysim,500,5,3,1,Nbigsim) b0mad = smad(ansmletr,dummysim,500,5,3,1,Nbigsim) print *, "The mean of the beta0 = ",b0mean print *, "The s.d. of the beta0 = ",b0sd print *, "The median of the beta0 = ",b0med print *, "The mad of the beta0 = ",b0mad print *, " " print *, "Results of the simulation, True MLE" print *, "beta1 = ",(ansmletr(i,4),i=1,Nbigsim) b1mean = smean(ansmletr,500,5,4,1,Nbigsim) b1sd = stmean(ansmletr,500,5,4,1,Nbigsim) b1med = smedian(ansmletr,dummysim,500,5,4,1,Nbigsim) b1mad = smad(ansmletr,dummysim,500,5,4,1,Nbigsim) b1mae = smae(ansmletr,dummysim,500,5,4,1,Nbigsim,0.83d0) b1mse = ((b1mean - 0.83) ** 2) + (b1sd **2) b1ase = ((b1med - 0.83) ** 2) + (b1mad **2) print *, "The mean of the beta1 = ",b1mean print *, "The s.d. of the beta1 = ",b1sd print *, "The mse of the beta1 = ",b1mse print *, "The mae of the beta1 = ",b1mae print *, "The median of the beta1 = ",b1med print *, "The mad of the beta1 = ",b1mad print *, "The madmse of the beta1 = ",b1ase print *, " " print *, "Results of the simulation, True MLE" print *, "sigma_e^2 = ",(ansmletr(i,1),i=1,Nbigsim) semean = smean(ansmletr,500,5,1,1,Nbigsim) sesd = stmean(ansmletr,500,5,1,1,Nbigsim) semed = smedian(ansmletr,dummysim,500,5,1,1,Nbigsim) semad = smad(ansmletr,dummysim,500,5,1,1,Nbigsim) print *, "The mean of the sigma_e^2 = ",semean print *, "The s.d. of the sigma_e^2 = ",sesd print *, "The median of the sigma_e^2 = ",semed print *, "The mad of the sigma_e^2 = ",semad print *, " " print *, "Results of the simulation, True MLE" print *, "sigma_u^2 = ",(ansmletr(i,2),i=1,Nbigsim) sumean = smean(ansmletr,500,5,2,1,Nbigsim) susd = stmean(ansmletr,500,5,2,1,Nbigsim) sumed = smedian(ansmletr,dummysim,500,5,2,1,Nbigsim) sumad = smad(ansmletr,dummysim,500,5,2,1,Nbigsim) print *, "The mean of the sigma_u^2 = ",sumean print *, "The s.d. of the sigma_u^2 = ",susd print *, "The median of the sigma_u^2 = ",sumed print *, "The mad of the sigma_u^2 = ",sumad print *, " " print *, "Results of the simulation, True MLE" print *, "rhoqt = ",(ansmletr(i,5),i=1,Nbigsim) rhomean = smean(ansmletr,500,5,5,1,Nbigsim) rhosd = stmean(ansmletr,500,5,5,1,Nbigsim) rhomed = smedian(ansmletr,dummysim,500,5,5,1,Nbigsim) rhomad = smad(ansmletr,dummysim,500,5,5,1,Nbigsim) rhomae = smae(ansmletr,dummysim,500,5,5,1,Nbigsim,0.54d0) rhomse = ((rhomean - 0.54) ** 2) + (rhosd **2) rhoase = ((rhomed - 0.54) ** 2) + (rhomad **2) print *, "The mean of the rho = ",rhomean print *, "The s.d. of the rho = ",rhosd print *, "The mse of the rho = ",rhomse print *, "The mae of the rho = ",rhomae print *, "The median of the rho = ",rhomed print *, "The mad of the rho = ",rhomad print *, "The madmse of the rho = ",rhoase do 455 i=1,500 cc = ansmletr(i,5) if (dabs(cc).ge.0.90d0) cc = 0.90 cc = dlog((1.0d0 + cc) / (1.0d0 - cc)) 455 ansmletr(i,5) = cc cc = dlog((1.0d0 + 0.54d0) / (1.0d0 - 0.54d0)) print *, " " print *, "Results of the simulation, True MLE" print *, "Fisher rhoqt = ",(ansmletr(i,5),i=1,Nbigsim) rrhomean = smean(ansmletr,500,5,5,1,Nbigsim) rrhosd = stmean(ansmletr,500,5,5,1,Nbigsim) rrhomed = smedian(ansmletr,dummysim,500,5,5,1,Nbigsim) rrhomad = smad(ansmletr,dummysim,500,5,5,1,Nbigsim) rrhomae = smae(ansmletr,dummysim,500,5,5,1,Nbigsim,cc) rrhomse = ((rrhomean - cc) ** 2) + (rrhosd **2) rrhoase = ((rrhomed - cc) ** 2) + (rrhomad **2) print *, "The mean of the Fisher rho = ",rrhomean print *, "The s.d. of the Fisher rho = ",rrhosd print *, "The mse of the Fisher rho = ",rrhomse print *, "The mae of the Fisher rho = ",rrhomae print *, "The median of the Fisher rho = ",rrhomed print *, "The mad of the Fisher rho = ",rrhomad print *, "The madmse of the Fisher rho = ",rrhoase smsumer(1,1) = b1mean smsumer(1,2) = b1sd smsumer(1,3) = b1med smsumer(1,4) = b1mad smsumer(1,5) = b1mae smsumer(1,6) = b1mse smsumer(1,7) = b1ase smsumer(2,1) = rhomean smsumer(2,2) = rhosd smsumer(2,3) = rhomed smsumer(2,4) = rhomad smsumer(2,5) = rhomae smsumer(2,6) = rhomse smsumer(2,7) = rhoase return end