2016-05-04 48 views
0

我正在研究具有不同参数(样本大小和方差)的正态分布和伽马分布的鲁棒性。模拟分布图

我也得到了模拟结果。他们是三张桌子。 但现在我不得不尝试绘制模拟分布图,以使人们更好地理解结果。

我在R中还是个新手。我是否需要在分布图中包含所有三个表的结果?

######################################################################## 
#For gamma distribution with equal skewness 1.5 

# rm(list=ls()) # clean the workspace 
nSims<-10000  #set the number of simulations 
alpha<-0.05  #set the significance level 

# to ensure the reproduction of the result 
# here we declare the random seed generator 
set.seed(1) 

#create vector to combine all std deviations 
sds<-matrix(c(4,4,6,4,8,4,10,4,12,4),nrow=2) 

sd1<-c(4,6,8,10,12) 
sd2<-c(4,4,4,4,4) 

## Put the samples sizes into matrix then use a loop for sample sizes 
sample_sizes<-matrix(c(10,10,10,25,25,25,25,50,25,100,50,25,50,100,100,25,100,100), 
nrow=2) 

#shape parameter for gamma distribution for equal skewness 
#forty five cases for each skewness!!!! 
sp1<-matrix(rep(c(16/9),each=45),ncol=1) 

scp <- c(1,1.5,2,2.5,3) 

##(use expand.grid)to create a data frame 
ss_scp<- expand.grid(sample_sizes[2,],scp) 

#create a matrix combining the forty five cases of combination of sample sizes,shape and scale parameter 
all <- cbind(rep(sample_sizes[1,], 5),ss_scp[,1],sp1,ss_scp[,2]) 

# name the column samples 1 and 2 and standard deviation 
colnames(all) <- c("m","n","sp","scp") 

#set empty vector of length no.of simulation(10000) to store p-value 
equal<-unequal<-mann<-c(rep(0,nrow(all))) 

#set nrow =nsims because wan storing every p-value simulated 
#for gamma distribution with equal skewness 
matrix_t <-matrix(0,nrow=nSims,ncol=5) 
matrix_u<-matrix(0,nrow=nSims,ncol=5) 
matrix_mann <-matrix(0,nrow=nSims,ncol=5) 

##for the samples sizes into matrix then use a loop for sample sizes 
# this loop steps through the all_combine matrix 
    for(ss in 1:nrow(all)) 
    { 
    #generate samples from the first column and second column 
    m<-all[ss,1] 
    n<-all[ss,2] 

     for (sim in 1:nSims) 
     { 
     #generate 2 random samples from gamma distribution with equal skewness 
     gamma1<-rgamma(m,all[ss,3],scale=all[ss,4]) 
     gamma2<-rgamma(n,all[ss,3],scale=1) 

     gamma1<-gamma1-all[ss,3]*all[ss,4] 
     gamma2<-gamma2-all[ss,3] 

     #extract p-value out and store every p-value into matrix 
     p<-t.test(gamma1,gamma2,var.equal=TRUE)$p.value 
     q<-t.test(gamma1,gamma2,var.equal=FALSE)$p.value 
     r<-wilcox.test(gamma1,gamma2)$p.value 

     matrix_t[sim,1]<- p 
     matrix_u[sim,1]<- q 
     matrix_mann[sim,1] <- r 
    } 
     ##store the result 
     equal[ss]<- sum(matrix_t[,1]<alpha) 
     unequal[ss]<-sum(matrix_u[,1]<alpha) 
     mann[ss]<- sum(matrix_mann[,1]<alpha) 
    } 

g1_equal<-cbind(all, equal, unequal, mann) 
print("g1_equal_skewness1.5)") 
print(g1_equal) 

    #samples sizes (10,10),(10,25).. 
    #standard deviation ratio (1,1.5,2,2.5,3) 
        Gamma(equal skewness)   Gamma(unequal skewness) 
        1.5 2.0 2.5 3.0 3.5  (1.5,1) 2,1.5 2.5,2 3,2.5 3.5,3 
10,10 
      Normal 
    1.0  506  382 379 343 270  246  422 426 383 303 247 
    1.5  472  493 463 507 537  571  531 518 548 528 532 
    2.0 516  597 679 736 829  935  597 680 760 836 951 
    2.5 498  627 747 905 1028 1215  687 825 944 1011 1197 
    3.0 493  678 864 1010 1190 1379  705 831 1015 1170 1436 


10,25 

    1.0 511  568 557 633 647 630  603 599 604  652 654 
    1.5 501  692 840 977 1012 1173  675 756 940  1068 1130 
    2.0 438  713 951 1049 1264 1470  773 869 1055 1259 1401 
    2.5 506  810 939 1101 1300 1594  761 960 1155 1339 1512 
    3.0 524  787 933 1176 1378 1599  772 967 1201 1339 1612 



25,25 

    1.0 479  463 451 447 417 414  513 429 439 469 392 
    1.5 493  534 556 504 568 587  537 517 528 539 555 
    2.0 510  543 599 676 663 773  538 607 677 712 725 
    2.5 487  591 662 731 807 908  581 643 733 769 893 
    3.0 488  614 668 761 811 1002 582 694 728 900 946 


25,50 

    1.0 519  585 487 569 559 579  521 572 568 581 583 
    1.5 510  532 651 695 725 836  625 647 729 737 802 
    2.0 501  586 660 758 846 888  618 653 794 876 957 
    2.5 466  635 687 823 937 996  612 702 782 909 1025 
    3.0 492  603 719 824 970 1045 640 704 826 945 1073 


25,100 

    1.0 486  559 589 670 726 778 552 614 666 752 750 
    1.5 494  621 700 787 903 955 602 703 774 842 1008 
    2.0 516  617 707 817 969 1073 613 755 774 932 1091 
    2.5 470  598 731 873 969 1118 624 752 849 970 1094 
    3.0 493  710 718 824 1021 1167 645 746 887 988 1149 

50,25 

    1.0 495  507 511 552 550 534  491 527 496 554 534 
    1.5 535  472 470 489 470 413  458 503 460 456 410 
    2.0 499  507 478 488 468 465  495 490 542 528 489 
    2.5 486  500 532 517 559 629  509 493 526 569 601 
    3.0 490  586 536 561 654 644  544 567 563 614 665 

50,100 1.0 518  515 530 531 514 569  516 494 517 548 578 
     1.5 528  503 542 597 596 656  554 565 612 606 708 
     2.0 453  525 588 640 727 775  520 625 628 727 772 
     2.5 500  586 660 669 733 837  552 622 660 695 802 
     3.0 494  557 640 680 747 847  582 634 686 776 834 

100,25             
     1.0 489  553 607 641 712 777  557 560 653 677 751 
     1.5 516  497 553 532 619 595  496 548 512 549 553 
     2.0 500  492 483 518 472 468  536 521 497 463 463 

     2.5 493  498 473 446 488 461  483 463 476 452 472 
     3.0 482  490 516 481 488 500  563 477 496 492 537 




100,100 
     1.0 472  508 492 483 517 487  517 521 476 505 485 
     1.5 507  498 496 511 518 546  520 520 498 547 531 
     2.0 465  478 540 542 584 599  496 504 585 558 589 
     2.5 508  486 566 551 614 602  520 539 583 601 642 
     3.0 494  497 575 545 614 651  561 557 590 615 624 
+1

SO帮助页面要求提供“最小”示例。这很难被称为“最小”,因为它需要很长时间来计算。目前我无法分辨您是否创建了无限循环或只是一个非常长时间运行的程序。我在6分钟后停止执行嵌套的for-loops。然后允许打印已经达到45行第34行的结果。你应该通过进一步编辑来改善你的问题,以解决这个问题,并且进一步的担忧是我不能说出你现在预期或暗示的那种“阴谋”。我只看到一张表,以便需要澄清。 –

+0

这只是其中的一种编码...我还有另一种编码...因此我有11个结果 – quess

+0

您似乎无法接受旨在提供帮助的建议。这是我努力将其锐化到一个很好的点。 1)缩短执行时间,然后2)说出需要的底模类型。 –

回答

0

这听起来像你想知道如何着手绘制你的分布。您可以创建三个独立的图,也可以将每个分布覆盖在一个图上。以下是如何在ggplot2中做到这一点。

假设您的数据框为df,列dist1,dist2dist3

install.packages('ggplot2') 
library(ggplot2) 

ggplot(df, aes(x = dist1)) + 
    geom_density(color = 'red') + 
    geom_density(aes(x = dist2), color = 'green') + 
    geom_density(aes(x = dist3), color = 'blue') 

这应该给你一个三行密度图,每个分布一个。如果你想创建三个独立的地块,只需为每个分配创建一个新的ggplot。

plot1 <- ggplot(df, aes(x = dist1) + geom_density() 
plot2 <- ggplot(df, aes(x = dist2) + geom_density() 

...等等。这有帮助吗?