2015年9月1日星期二
2015年8月16日星期日
Varince-covariance calcuation for logistic regression
library(Matrix)
library(sandwich)
mydata <- read.csv(http://www.ats.ucla.edu/stat/data/binary.csv)
mylogit <- glm(admit ~ gre + gpa, data = mydata, family = "binomial")
X <- as.matrix(cbind(1, mydata[,c('gre','gpa')]))
n <- nrow(X)
pi<-mylogit$fit
w<-pi*(1-pi)
v<-Diagonal(n, x = w)
var_b<-solve(t(X)%*%v%*%X)
var_b
vcov(mylogit)
library(sandwich)
mydata <- read.csv(http://www.ats.ucla.edu/stat/data/binary.csv)
mylogit <- glm(admit ~ gre + gpa, data = mydata, family = "binomial")
X <- as.matrix(cbind(1, mydata[,c('gre','gpa')]))
n <- nrow(X)
pi<-mylogit$fit
w<-pi*(1-pi)
v<-Diagonal(n, x = w)
var_b<-solve(t(X)%*%v%*%X)
var_b
vcov(mylogit)
2015年8月7日星期五
Understanding and correcting complete or quasi-complete separation problems
This is a common problem with logistic models. Unlike in ordinary least-squares regression for modeling a normally distributed response, when a logistic model perfectly or nearly perfectly predicts the response (that is, separates the response levels), unique maximum likelihood estimates do not exist. Some model parameters are actually infinite. This is a common result of the data being sparse, meaning that not all response levels are observed in each of the predictor settings, which often happens with small data sets or when the event is rare. See "Existence of Maximum Likelihood Estimates" in the Details section of the LOGISTIC documentation.
When the separation condition is detected, the procedure terminates the maximum likelihood iteration process and reports the results at the last iteration. The resulting model is often quite good at classifying observations and might be usable for this purpose, but inferences about the parameters based on Wald chi-square statistics should be avoided. For details about the stopping rule, see the "Existence of Maximum Likelihood Estimates" section of the LOGISTIC documentation and the So reference.
One possible solution to the separation problem is to make the model fit less perfectly by reducing the number of variables or effects in the model, categorizing continuous variables, or merging categories of categorical (CLASS) variables. It is often difficult to know exactly which variables cause the separation, but variables that exhibit large parameter estimates or standard errors are likely candidates.
When developing a model, it can be very helpful to use forward or stepwise (not backward) selection via the SELECTION= option in the MODEL statement. These methods build a model that might avoid separation, or if separation occurs in a model fit at a particular step, you might consider using the model successfully fit at the previous step.
Valid tests of model effects can be obtained by using likelihood ratio (LR) tests rather than the Wald tests that are presented by PROC LOGISTIC. For models that do not involve survey data, the easiest way to do this is to fit the model in PROC GENMOD (use the DIST=BINOMIAL option in the MODEL statement) and use the TYPE3 option or the CONTRAST statement to obtain LR tests.
Because separation often happens with small data sets, another possible solution is to use exact logistic regression, which is available via the EXACT statement in PROC LOGISTIC. However, exact methods are very computationally intensive and can take great amounts of time and/or memory even for moderate sized problems. A penalized likelihood method by Firth can often converge to a solution even though a unique solution by maximum likelihood does not exist in separation cases. Firth's method is available, beginning in SAS 9.2, by specifying the FIRTH option in the MODEL statement.
If interest lies exclusively in testing the association of one predictor with the response while adjusting for the effects of the other predictors, then you could use the CMH option in PROC FREQ. All predictors other than the one of interest are used as stratifying variables. Since no model is fit by PROC FREQ, the problem of infinite parameters caused by separation is not an issue. For example, consider the final model involving Sex, Age, and Treatment in the example titled "Logistic Modeling with Categorical Predictors" in the LOGISTIC documentation. The following statements test the association of Treatment with the response, Pain, while adjusting for Sex and Age. The NOPRINT option is used to avoid displaying the large multi-way table. The CMH option computes the test. Since Treatment is an unordered categorical variable with three levels, the second CMH statistic labeled "Row Mean Scores Differ" is an appropriate test.
Usage Note 22599: Understanding and correcting complete or quasi-complete separation problems
When the separation condition is detected, the procedure terminates the maximum likelihood iteration process and reports the results at the last iteration. The resulting model is often quite good at classifying observations and might be usable for this purpose, but inferences about the parameters based on Wald chi-square statistics should be avoided. For details about the stopping rule, see the "Existence of Maximum Likelihood Estimates" section of the LOGISTIC documentation and the So reference.
One possible solution to the separation problem is to make the model fit less perfectly by reducing the number of variables or effects in the model, categorizing continuous variables, or merging categories of categorical (CLASS) variables. It is often difficult to know exactly which variables cause the separation, but variables that exhibit large parameter estimates or standard errors are likely candidates.
When developing a model, it can be very helpful to use forward or stepwise (not backward) selection via the SELECTION= option in the MODEL statement. These methods build a model that might avoid separation, or if separation occurs in a model fit at a particular step, you might consider using the model successfully fit at the previous step.
Valid tests of model effects can be obtained by using likelihood ratio (LR) tests rather than the Wald tests that are presented by PROC LOGISTIC. For models that do not involve survey data, the easiest way to do this is to fit the model in PROC GENMOD (use the DIST=BINOMIAL option in the MODEL statement) and use the TYPE3 option or the CONTRAST statement to obtain LR tests.
Because separation often happens with small data sets, another possible solution is to use exact logistic regression, which is available via the EXACT statement in PROC LOGISTIC. However, exact methods are very computationally intensive and can take great amounts of time and/or memory even for moderate sized problems. A penalized likelihood method by Firth can often converge to a solution even though a unique solution by maximum likelihood does not exist in separation cases. Firth's method is available, beginning in SAS 9.2, by specifying the FIRTH option in the MODEL statement.
If interest lies exclusively in testing the association of one predictor with the response while adjusting for the effects of the other predictors, then you could use the CMH option in PROC FREQ. All predictors other than the one of interest are used as stratifying variables. Since no model is fit by PROC FREQ, the problem of infinite parameters caused by separation is not an issue. For example, consider the final model involving Sex, Age, and Treatment in the example titled "Logistic Modeling with Categorical Predictors" in the LOGISTIC documentation. The following statements test the association of Treatment with the response, Pain, while adjusting for Sex and Age. The NOPRINT option is used to avoid displaying the large multi-way table. The CMH option computes the test. Since Treatment is an unordered categorical variable with three levels, the second CMH statistic labeled "Row Mean Scores Differ" is an appropriate test.
proc freq data=Neuralgia;
table Sex * Age * Treatment * Pain / noprint cmh;
run;
The resulting p-value (0.0040) closely matches the p-value of the Treatment Type3 test from the logistic model (0.0018) Usage Note 22599: Understanding and correcting complete or quasi-complete separation problems
2015年6月30日星期二
2015年6月27日星期六
EM alogrithm for mixed Normal distribtuion
Data are from Robert Hogg Exercise 6.6.8
library(car)
x<-c(119.0,96.0,146.2,138.6,143.4,98.2,124.5,114.1,136.2,136.4,
184.8,79.8,151.9,114.2,145.7,95.9,97.3,136.4,109.2,103.2)
densityPlot(x)
p_0<-0.6
u1_0<-105
u2_0<-130
sigma1_0<-15
sigma2_0<-25
p<-c()
pi<-c()
u1<-c()
u2<-c()
sigma1_2<-c()
sigma2_2<-c()
s<-c()
s1<-c()
diff<-c()
j<-1 #this is imporant
s<-c(p_0,u1_0,u2_0,sigma1_0^2,sigma2_0^2) #initial guess of p,u1,u2,sigma1, sigma2
p_i_0<-p_0*dnorm(x,u2_0,sigma2_0)/((1-p_0)*dnorm(x,u1_0,sigma1_0)+p_0*dnorm(x,u2_0,sigma2_0))
p[j]<-mean(p_i_0)
u1[j]<-sum((1-p_i_0)*x)/sum(1-p_i_0)
u2[j]<-sum(p_i_0*x)/sum(p_i_0)
sigma1_2[j]<-sum((1-p_i_0)*(x-u1_0)^2)/sum(1-p_i_0)
sigma2_2[j]<-sum(p_i_0*(x-u2_0)^2)/sum(p_i_0)
s1<-c(p[j],u1[j],u2[j],sigma1_2[j],sigma2_2[j])
for(i in 1:5){
diff[i]<-abs(s1[i]-s[i])
}
for(i in 1:5){ #in fact we don't need this, we can mainly focus on p
while(diff[i]>0.000001){
s<-c(p[j],u1[j],u2[j],sigma1_2[j],sigma2_2[j])
p_i<-p[j]*dnorm(x,u2[j],(sigma2_2[j])^0.5)/((1-p[j])*dnorm(x,u1[j],sigma1_2[j]^0.5)+p[j]*dnorm(x,u2[j],sigma2_2[j]^0.5))
#p_i will be updated by p,u1,u2,sigma1 and sigma2
u1[j+1]<-sum((1-p_i)*x)/sum(1-p_i)
u2[j+1]<-sum(p_i*x)/sum(p_i)
sigma1_2[j+1]<-sum((1-p_i)*(x-u1[j+1])^2)/sum(1-p_i)
sigma2_2[j+1]<-sum(p_i*(x-u2[j+1])^2)/sum(p_i)
p[j+1]<-mean(p_i)
s1<-c(p[j+1],u1[j+1],u2[j+1],sigma1_2[j+1],sigma2_2[j+1])
diff[i]<-abs(s1[i]-s[i])
j=j+1
}
}
#results:
After 301 iterations all parameters are within 0.00000 accuracy
For p:
[1] 0.6446007 0.6508217 0.6549901 0.6582507 0.6608042 0.6629207 0.6649003
[8] 0.6668472 0.6687443 0.6705665 0.6723035 0.6739539 0.6755194 0.6770028
[15] 0.6784075 0.6797370 0.6809951 0.6821853 0.6833114 0.6843770 0.6853854
[22] 0.6863398 0.6872434 0.6880990 0.6889093 0.6896769 0.6904043 0.6910936
[29] 0.6917470 0.6923664 0.6929539 0.6935110 0.6940395 0.6945410 0.6950169
[36] 0.6954685 0.6958972 0.6963042 0.6966907 0.6970577 0.6974063 0.6977374
[43] 0.6980519 0.6983508 0.6986348 0.6989047 0.6991612 0.6994050 0.6996367
[50] 0.6998571 0.7000666 0.7002657 0.7004552 0.7006353 0.7008066 0.7009696
[57] 0.7011245 0.7012720 0.7014122 0.7015457 0.7016726 0.7017934 0.7019083
[64] 0.7020177 0.7021217 0.7022208 0.7023150 0.7024047 0.7024901 0.7025713
[71] 0.7026487 0.7027223 0.7027924 0.7028591 0.7029226 0.7029830 0.7030405
[78] 0.7030953 0.7031475 0.7031971 0.7032444 0.7032894 0.7033322 0.7033730
[85] 0.7034119 0.7034489 0.7034841 0.7035176 0.7035496 0.7035800 0.7036089
[92] 0.7036365 0.7036627 0.7036877 0.7037115 0.7037342 0.7037558 0.7037764
[99] 0.7037959 0.7038146 0.7038323 0.7038493 0.7038654 0.7038807 0.7038953
[106] 0.7039092 0.7039225 0.7039351 0.7039471 0.7039585 0.7039694 0.7039798
[113] 0.7039897 0.7039991 0.7040081 0.7040166 0.7040247 0.7040325 0.7040399
[120] 0.7040469 0.7040536 0.7040600 0.7040660 0.7040718 0.7040773 0.7040826
[127] 0.7040876 0.7040923 0.7040968 0.7041012 0.7041053 0.7041092 0.7041129
[134] 0.7041165 0.7041199 0.7041231 0.7041261 0.7041291 0.7041319 0.7041345
[141] 0.7041370 0.7041394 0.7041417 0.7041439 0.7041460 0.7041480 0.7041499
[148] 0.7041516 0.7041534 0.7041550 0.7041565 0.7041580 0.7041594 0.7041608
[155] 0.7041620 0.7041633 0.7041644 0.7041655 0.7041666 0.7041676 0.7041685
[162] 0.7041694 0.7041703 0.7041711 0.7041719 0.7041727 0.7041734 0.7041740
[169] 0.7041747 0.7041753 0.7041759 0.7041765 0.7041770 0.7041775 0.7041780
[176] 0.7041784 0.7041789 0.7041793 0.7041797 0.7041801 0.7041804 0.7041808
[183] 0.7041811 0.7041814 0.7041817 0.7041820 0.7041822 0.7041825 0.7041827
[190] 0.7041830 0.7041832 0.7041834 0.7041836 0.7041838 0.7041840 0.7041842
[197] 0.7041843 0.7041845 0.7041846 0.7041848 0.7041849 0.7041850 0.7041852
[204] 0.7041853 0.7041854 0.7041855 0.7041856 0.7041857 0.7041858 0.7041859
[211] 0.7041860 0.7041860 0.7041861 0.7041862 0.7041863 0.7041863 0.7041864
[218] 0.7041864 0.7041865 0.7041866 0.7041866 0.7041867 0.7041867 0.7041867
[225] 0.7041868 0.7041868 0.7041869 0.7041869 0.7041869 0.7041870 0.7041870
[232] 0.7041870 0.7041871 0.7041871 0.7041871 0.7041871 0.7041872 0.7041872
[239] 0.7041872 0.7041872 0.7041872 0.7041873 0.7041873 0.7041873 0.7041873
[246] 0.7041873 0.7041873 0.7041874 0.7041874 0.7041874 0.7041874 0.7041874
[253] 0.7041874 0.7041874 0.7041874 0.7041875 0.7041875 0.7041875 0.7041875
[260] 0.7041875 0.7041875 0.7041875 0.7041875 0.7041875 0.7041875 0.7041875
[267] 0.7041875 0.7041875 0.7041875 0.7041875 0.7041875 0.7041876 0.7041876
[274] 0.7041876 0.7041876 0.7041876 0.7041876 0.7041876 0.7041876 0.7041876
[281] 0.7041876 0.7041876 0.7041876 0.7041876 0.7041876 0.7041876 0.7041876
[288] 0.7041876 0.7041876 0.7041876 0.7041876 0.7041876 0.7041876 0.7041876
[295] 0.7041876 0.7041876 0.7041876 0.7041876 0.7041876 0.7041876 0.7041876
For u1
u1
[1] 105.37122 104.28497 103.19437 102.19782 101.34808 100.69007 100.21667 99.88967 99.66755
[10] 99.51631 99.41138 99.33623 99.28022 99.23661 99.20121 99.17139 99.14549 99.12247
[19] 99.10163 99.08255 99.06490 99.04849 99.03315 99.01877 99.00527 98.99257 98.98060
[28] 98.96932 98.95867 98.94862 98.93912 98.93015 98.92166 98.91364 98.90605 98.89886
[37] 98.89206 98.88562 98.87952 98.87374 98.86826 98.86307 98.85815 98.85348 98.84905
[46] 98.84485 98.84087 98.83709 98.83350 98.83009 98.82686 98.82379 98.82087 98.81810
[55] 98.81547 98.81297 98.81059 98.80833 98.80619 98.80415 98.80221 98.80037 98.79861
[64] 98.79695 98.79536 98.79386 98.79243 98.79106 98.78977 98.78853 98.78736 98.78624
[73] 98.78518 98.78417 98.78321 98.78230 98.78143 98.78060 98.77981 98.77906 98.77835
[82] 98.77767 98.77702 98.77640 98.77582 98.77526 98.77473 98.77422 98.77374 98.77329
[91] 98.77285 98.77243 98.77204 98.77166 98.77130 98.77096 98.77064 98.77033 98.77003
[100] 98.76975 98.76949 98.76923 98.76899 98.76876 98.76854 98.76833 98.76813 98.76794
[109] 98.76776 98.76759 98.76743 98.76727 98.76712 98.76698 98.76685 98.76672 98.76660
[118] 98.76648 98.76637 98.76626 98.76616 98.76607 98.76598 98.76589 98.76581 98.76573
[127] 98.76565 98.76558 98.76551 98.76545 98.76539 98.76533 98.76527 98.76522 98.76517
[136] 98.76512 98.76507 98.76503 98.76499 98.76495 98.76491 98.76487 98.76484 98.76481
[145] 98.76478 98.76475 98.76472 98.76469 98.76467 98.76464 98.76462 98.76460 98.76457
[154] 98.76455 98.76454 98.76452 98.76450 98.76448 98.76447 98.76445 98.76444 98.76442
[163] 98.76441 98.76440 98.76439 98.76438 98.76437 98.76436 98.76435 98.76434 98.76433
[172] 98.76432 98.76431 98.76430 98.76430 98.76429 98.76428 98.76428 98.76427 98.76427
[181] 98.76426 98.76425 98.76425 98.76425 98.76424 98.76424 98.76423 98.76423 98.76423
[190] 98.76422 98.76422 98.76422 98.76421 98.76421 98.76421 98.76420 98.76420 98.76420
[199] 98.76420 98.76419 98.76419 98.76419 98.76419 98.76419 98.76419 98.76418 98.76418
[208] 98.76418 98.76418 98.76418 98.76418 98.76418 98.76417 98.76417 98.76417 98.76417
[217] 98.76417 98.76417 98.76417 98.76417 98.76417 98.76417 98.76417 98.76417 98.76416
[226] 98.76416 98.76416 98.76416 98.76416 98.76416 98.76416 98.76416 98.76416 98.76416
[235] 98.76416 98.76416 98.76416 98.76416 98.76416 98.76416 98.76416 98.76416 98.76416
[244] 98.76416 98.76416 98.76416 98.76416 98.76416 98.76416 98.76416 98.76416 98.76416
[253] 98.76415 98.76415 98.76415 98.76415 98.76415 98.76415 98.76415 98.76415 98.76415
[262] 98.76415 98.76415 98.76415 98.76415 98.76415 98.76415 98.76415 98.76415 98.76415
[271] 98.76415 98.76415 98.76415 98.76415 98.76415 98.76415 98.76415 98.76415 98.76415
[280] 98.76415 98.76415 98.76415 98.76415 98.76415 98.76415 98.76415 98.76415 98.76415
[289] 98.76415 98.76415 98.76415 98.76415 98.76415 98.76415 98.76415 98.76415 98.76415
[298] 98.76415 98.76415 98.76415 98.76415
for u2
[1] 133.5728 133.8861 134.2721 134.6356 134.9464 135.1737 135.3096 135.3706 135.3799 135.3572
[11] 135.3157 135.2642 135.2078 135.1499 135.0923 135.0362 134.9820 134.9302 134.8809 134.8340
[21] 134.7895 134.7473 134.7073 134.6695 134.6337 134.5997 134.5676 134.5371 134.5083 134.4809
[31] 134.4550 134.4305 134.4072 134.3851 134.3641 134.3442 134.3254 134.3075 134.2905 134.2743
[41] 134.2590 134.2445 134.2307 134.2175 134.2051 134.1932 134.1820 134.1713 134.1611 134.1515
[51] 134.1423 134.1335 134.1252 134.1173 134.1098 134.1027 134.0959 134.0895 134.0833 134.0775
[61] 134.0719 134.0666 134.0616 134.0568 134.0523 134.0479 134.0438 134.0399 134.0361 134.0326
[71] 134.0292 134.0260 134.0229 134.0200 134.0172 134.0146 134.0121 134.0097 134.0074 134.0052
[81] 134.0032 134.0012 133.9993 133.9975 133.9958 133.9942 133.9927 133.9912 133.9898 133.9885
[91] 133.9872 133.9860 133.9849 133.9838 133.9827 133.9818 133.9808 133.9799 133.9791 133.9782
[101] 133.9775 133.9767 133.9760 133.9754 133.9747 133.9741 133.9735 133.9730 133.9725 133.9720
[111] 133.9715 133.9710 133.9706 133.9702 133.9698 133.9694 133.9691 133.9687 133.9684 133.9681
[121] 133.9678 133.9675 133.9673 133.9670 133.9668 133.9665 133.9663 133.9661 133.9659 133.9657
[131] 133.9655 133.9654 133.9652 133.9651 133.9649 133.9648 133.9646 133.9645 133.9644 133.9643
[141] 133.9642 133.9640 133.9639 133.9639 133.9638 133.9637 133.9636 133.9635 133.9634 133.9634
[151] 133.9633 133.9632 133.9632 133.9631 133.9631 133.9630 133.9630 133.9629 133.9629 133.9628
[161] 133.9628 133.9627 133.9627 133.9627 133.9626 133.9626 133.9626 133.9625 133.9625 133.9625
[171] 133.9625 133.9624 133.9624 133.9624 133.9624 133.9623 133.9623 133.9623 133.9623 133.9623
[181] 133.9623 133.9622 133.9622 133.9622 133.9622 133.9622 133.9622 133.9622 133.9622 133.9621
[191] 133.9621 133.9621 133.9621 133.9621 133.9621 133.9621 133.9621 133.9621 133.9621 133.9621
[201] 133.9621 133.9621 133.9620 133.9620 133.9620 133.9620 133.9620 133.9620 133.9620 133.9620
[211] 133.9620 133.9620 133.9620 133.9620 133.9620 133.9620 133.9620 133.9620 133.9620 133.9620
[221] 133.9620 133.9620 133.9620 133.9620 133.9620 133.9620 133.9620 133.9620 133.9620 133.9620
[231] 133.9620 133.9620 133.9620 133.9620 133.9620 133.9620 133.9620 133.9620 133.9620 133.9620
[241] 133.9620 133.9620 133.9620 133.9620 133.9620 133.9620 133.9620 133.9620 133.9620 133.9620
[251] 133.9620 133.9620 133.9620 133.9620 133.9619 133.9619 133.9619 133.9619 133.9619 133.9619
[261] 133.9619 133.9619 133.9619 133.9619 133.9619 133.9619 133.9619 133.9619 133.9619 133.9619
[271] 133.9619 133.9619 133.9619 133.9619 133.9619 133.9619 133.9619 133.9619 133.9619 133.9619
[281] 133.9619 133.9619 133.9619 133.9619 133.9619 133.9619 133.9619 133.9619 133.9619 133.9619
[291] 133.9619 133.9619 133.9619 133.9619 133.9619 133.9619 133.9619 133.9619 133.9619 133.9619
[301] 133.9619
For sigma1 square:
[1] 220.12489 205.10096 183.53226 160.65836 141.18629 127.45956 118.93847 113.86190 110.77411
[10] 108.80379 107.47407 106.52342 105.80421 105.23069 104.75194 104.33723 103.96770 103.63162
[19] 103.32155 103.03264 102.76167 102.50639 102.26517 102.03679 101.82027 101.61481 101.41971
[28] 101.23438 101.05826 100.89085 100.73171 100.58039 100.43649 100.29964 100.16948 100.04567
[37] 99.92790 99.81586 99.70927 99.60785 99.51136 99.41954 99.33217 99.24902 99.16990
[46] 99.09459 99.02292 98.95471 98.88979 98.82799 98.76917 98.71318 98.65988 98.60914
[55] 98.56084 98.51486 98.47108 98.42941 98.38973 98.35196 98.31599 98.28174 98.24914
[64] 98.21810 98.18854 98.16039 98.13359 98.10807 98.08377 98.06063 98.03860 98.01762
[73] 97.99764 97.97861 97.96049 97.94324 97.92681 97.91116 97.89626 97.88207 97.86855
[82] 97.85568 97.84342 97.83175 97.82063 97.81005 97.79996 97.79036 97.78122 97.77251
[91] 97.76421 97.75631 97.74879 97.74162 97.73479 97.72829 97.72210 97.71621 97.71059
[100] 97.70524 97.70015 97.69530 97.69067 97.68627 97.68208 97.67809 97.67429 97.67067
[109] 97.66722 97.66393 97.66080 97.65782 97.65498 97.65228 97.64970 97.64725 97.64491
[118] 97.64269 97.64057 97.63855 97.63663 97.63480 97.63305 97.63139 97.62981 97.62830
[127] 97.62687 97.62550 97.62420 97.62296 97.62178 97.62065 97.61958 97.61856 97.61759
[136] 97.61666 97.61578 97.61494 97.61414 97.61338 97.61265 97.61196 97.61130 97.61067
[145] 97.61008 97.60951 97.60897 97.60845 97.60796 97.60749 97.60704 97.60662 97.60621
[154] 97.60583 97.60546 97.60511 97.60478 97.60446 97.60416 97.60387 97.60360 97.60334
[163] 97.60309 97.60285 97.60263 97.60241 97.60221 97.60201 97.60183 97.60165 97.60148
[172] 97.60132 97.60117 97.60102 97.60088 97.60075 97.60062 97.60051 97.60039 97.60028
[181] 97.60018 97.60008 97.59999 97.59990 97.59981 97.59973 97.59965 97.59958 97.59951
[190] 97.59944 97.59938 97.59932 97.59926 97.59921 97.59915 97.59910 97.59906 97.59901
[199] 97.59897 97.59893 97.59889 97.59885 97.59881 97.59878 97.59875 97.59872 97.59869
[208] 97.59866 97.59863 97.59861 97.59859 97.59856 97.59854 97.59852 97.59850 97.59848
[217] 97.59846 97.59845 97.59843 97.59841 97.59840 97.59839 97.59837 97.59836 97.59835
[226] 97.59834 97.59832 97.59831 97.59830 97.59829 97.59829 97.59828 97.59827 97.59826
[235] 97.59825 97.59825 97.59824 97.59823 97.59823 97.59822 97.59822 97.59821 97.59821
[244] 97.59820 97.59820 97.59819 97.59819 97.59818 97.59818 97.59818 97.59817 97.59817
[253] 97.59817 97.59816 97.59816 97.59816 97.59816 97.59815 97.59815 97.59815 97.59815
[262] 97.59814 97.59814 97.59814 97.59814 97.59814 97.59814 97.59813 97.59813 97.59813
[271] 97.59813 97.59813 97.59813 97.59813 97.59813 97.59812 97.59812 97.59812 97.59812
[280] 97.59812 97.59812 97.59812 97.59812 97.59812 97.59812 97.59812 97.59812 97.59812
[289] 97.59811 97.59811 97.59811 97.59811 97.59811 97.59811 97.59811 97.59811 97.59811
[298] 97.59811 97.59811 97.59811 97.59811
For sigma2 square
[1] 559.5878 525.6849 505.7973 488.0545 472.0863 458.8586 449.1209 442.7654 439.1009 437.3237
[11] 436.7687 436.9631 437.5961 438.4720 439.4702 440.5184 441.5732 442.6096 443.6137 444.5780
[21] 445.4991 446.3760 447.2092 447.9998 448.7495 449.4602 450.1336 450.7719 451.3767 451.9501
[31] 452.4935 453.0088 453.4975 453.9609 454.4006 454.8177 455.2135 455.5891 455.9457 456.2843
[41] 456.6057 456.9110 457.2010 457.4764 457.7380 457.9866 458.2229 458.4474 458.6608 458.8636
[51] 459.0564 459.2397 459.4140 459.5797 459.7373 459.8872 460.0297 460.1653 460.2943 460.4170
[61] 460.5337 460.6447 460.7503 460.8509 460.9465 461.0375 461.1241 461.2065 461.2850 461.3596
[71] 461.4307 461.4983 461.5627 461.6240 461.6823 461.7378 461.7907 461.8410 461.8889 461.9345
[81] 461.9779 462.0192 462.0586 462.0960 462.1317 462.1657 462.1980 462.2288 462.2581 462.2861
[91] 462.3126 462.3380 462.3621 462.3850 462.4069 462.4277 462.4475 462.4664 462.4844 462.5015
[101] 462.5178 462.5333 462.5481 462.5622 462.5756 462.5883 462.6005 462.6121 462.6231 462.6336
[111] 462.6436 462.6532 462.6622 462.6709 462.6791 462.6869 462.6944 462.7015 462.7083 462.7147
[121] 462.7209 462.7267 462.7323 462.7376 462.7427 462.7475 462.7521 462.7564 462.7606 462.7646
[131] 462.7683 462.7719 462.7754 462.7786 462.7817 462.7847 462.7875 462.7902 462.7927 462.7952
[141] 462.7975 462.7997 462.8018 462.8038 462.8057 462.8075 462.8093 462.8109 462.8125 462.8140
[151] 462.8154 462.8167 462.8180 462.8193 462.8204 462.8216 462.8226 462.8236 462.8246 462.8255
[161] 462.8264 462.8272 462.8280 462.8288 462.8295 462.8302 462.8308 462.8315 462.8321 462.8326
[171] 462.8332 462.8337 462.8342 462.8346 462.8351 462.8355 462.8359 462.8363 462.8366 462.8370
[181] 462.8373 462.8376 462.8379 462.8382 462.8385 462.8387 462.8390 462.8392 462.8394 462.8397
[191] 462.8399 462.8401 462.8402 462.8404 462.8406 462.8407 462.8409 462.8410 462.8412 462.8413
[201] 462.8414 462.8416 462.8417 462.8418 462.8419 462.8420 462.8421 462.8422 462.8422 462.8423
[211] 462.8424 462.8425 462.8425 462.8426 462.8427 462.8427 462.8428 462.8428 462.8429 462.8429
[221] 462.8430 462.8430 462.8431 462.8431 462.8432 462.8432 462.8432 462.8433 462.8433 462.8433
[231] 462.8434 462.8434 462.8434 462.8434 462.8435 462.8435 462.8435 462.8435 462.8435 462.8436
[241] 462.8436 462.8436 462.8436 462.8436 462.8436 462.8437 462.8437 462.8437 462.8437 462.8437
[251] 462.8437 462.8437 462.8437 462.8437 462.8438 462.8438 462.8438 462.8438 462.8438 462.8438
[261] 462.8438 462.8438 462.8438 462.8438 462.8438 462.8438 462.8438 462.8438 462.8438 462.8438
[271] 462.8439 462.8439 462.8439 462.8439 462.8439 462.8439 462.8439 462.8439 462.8439 462.8439
[281] 462.8439 462.8439 462.8439 462.8439 462.8439 462.8439 462.8439 462.8439 462.8439 462.8439
[291] 462.8439 462.8439 462.8439 462.8439 462.8439 462.8439 462.8439 462.8439 462.8439 462.8439
[301] 462.8439
library(car)
x<-c(119.0,96.0,146.2,138.6,143.4,98.2,124.5,114.1,136.2,136.4,
184.8,79.8,151.9,114.2,145.7,95.9,97.3,136.4,109.2,103.2)
densityPlot(x)
p_0<-0.6
u1_0<-105
u2_0<-130
sigma1_0<-15
sigma2_0<-25
p<-c()
pi<-c()
u1<-c()
u2<-c()
sigma1_2<-c()
sigma2_2<-c()
s<-c()
s1<-c()
diff<-c()
j<-1 #this is imporant
s<-c(p_0,u1_0,u2_0,sigma1_0^2,sigma2_0^2) #initial guess of p,u1,u2,sigma1, sigma2
p_i_0<-p_0*dnorm(x,u2_0,sigma2_0)/((1-p_0)*dnorm(x,u1_0,sigma1_0)+p_0*dnorm(x,u2_0,sigma2_0))
p[j]<-mean(p_i_0)
u1[j]<-sum((1-p_i_0)*x)/sum(1-p_i_0)
u2[j]<-sum(p_i_0*x)/sum(p_i_0)
sigma1_2[j]<-sum((1-p_i_0)*(x-u1_0)^2)/sum(1-p_i_0)
sigma2_2[j]<-sum(p_i_0*(x-u2_0)^2)/sum(p_i_0)
s1<-c(p[j],u1[j],u2[j],sigma1_2[j],sigma2_2[j])
for(i in 1:5){
diff[i]<-abs(s1[i]-s[i])
}
for(i in 1:5){ #in fact we don't need this, we can mainly focus on p
while(diff[i]>0.000001){
s<-c(p[j],u1[j],u2[j],sigma1_2[j],sigma2_2[j])
p_i<-p[j]*dnorm(x,u2[j],(sigma2_2[j])^0.5)/((1-p[j])*dnorm(x,u1[j],sigma1_2[j]^0.5)+p[j]*dnorm(x,u2[j],sigma2_2[j]^0.5))
#p_i will be updated by p,u1,u2,sigma1 and sigma2
u1[j+1]<-sum((1-p_i)*x)/sum(1-p_i)
u2[j+1]<-sum(p_i*x)/sum(p_i)
sigma1_2[j+1]<-sum((1-p_i)*(x-u1[j+1])^2)/sum(1-p_i)
sigma2_2[j+1]<-sum(p_i*(x-u2[j+1])^2)/sum(p_i)
p[j+1]<-mean(p_i)
s1<-c(p[j+1],u1[j+1],u2[j+1],sigma1_2[j+1],sigma2_2[j+1])
diff[i]<-abs(s1[i]-s[i])
j=j+1
}
}
#results:
For p:
[1] 0.6446007 0.6508217 0.6549901 0.6582507 0.6608042 0.6629207 0.6649003
[8] 0.6668472 0.6687443 0.6705665 0.6723035 0.6739539 0.6755194 0.6770028
[15] 0.6784075 0.6797370 0.6809951 0.6821853 0.6833114 0.6843770 0.6853854
[22] 0.6863398 0.6872434 0.6880990 0.6889093 0.6896769 0.6904043 0.6910936
[29] 0.6917470 0.6923664 0.6929539 0.6935110 0.6940395 0.6945410 0.6950169
[36] 0.6954685 0.6958972 0.6963042 0.6966907 0.6970577 0.6974063 0.6977374
[43] 0.6980519 0.6983508 0.6986348 0.6989047 0.6991612 0.6994050 0.6996367
[50] 0.6998571 0.7000666 0.7002657 0.7004552 0.7006353 0.7008066 0.7009696
[57] 0.7011245 0.7012720 0.7014122 0.7015457 0.7016726 0.7017934 0.7019083
[64] 0.7020177 0.7021217 0.7022208 0.7023150 0.7024047 0.7024901 0.7025713
[71] 0.7026487 0.7027223 0.7027924 0.7028591 0.7029226 0.7029830 0.7030405
[78] 0.7030953 0.7031475 0.7031971 0.7032444 0.7032894 0.7033322 0.7033730
[85] 0.7034119 0.7034489 0.7034841 0.7035176 0.7035496 0.7035800 0.7036089
[92] 0.7036365 0.7036627 0.7036877 0.7037115 0.7037342 0.7037558 0.7037764
[99] 0.7037959 0.7038146 0.7038323 0.7038493 0.7038654 0.7038807 0.7038953
[106] 0.7039092 0.7039225 0.7039351 0.7039471 0.7039585 0.7039694 0.7039798
[113] 0.7039897 0.7039991 0.7040081 0.7040166 0.7040247 0.7040325 0.7040399
[120] 0.7040469 0.7040536 0.7040600 0.7040660 0.7040718 0.7040773 0.7040826
[127] 0.7040876 0.7040923 0.7040968 0.7041012 0.7041053 0.7041092 0.7041129
[134] 0.7041165 0.7041199 0.7041231 0.7041261 0.7041291 0.7041319 0.7041345
[141] 0.7041370 0.7041394 0.7041417 0.7041439 0.7041460 0.7041480 0.7041499
[148] 0.7041516 0.7041534 0.7041550 0.7041565 0.7041580 0.7041594 0.7041608
[155] 0.7041620 0.7041633 0.7041644 0.7041655 0.7041666 0.7041676 0.7041685
[162] 0.7041694 0.7041703 0.7041711 0.7041719 0.7041727 0.7041734 0.7041740
[169] 0.7041747 0.7041753 0.7041759 0.7041765 0.7041770 0.7041775 0.7041780
[176] 0.7041784 0.7041789 0.7041793 0.7041797 0.7041801 0.7041804 0.7041808
[183] 0.7041811 0.7041814 0.7041817 0.7041820 0.7041822 0.7041825 0.7041827
[190] 0.7041830 0.7041832 0.7041834 0.7041836 0.7041838 0.7041840 0.7041842
[197] 0.7041843 0.7041845 0.7041846 0.7041848 0.7041849 0.7041850 0.7041852
[204] 0.7041853 0.7041854 0.7041855 0.7041856 0.7041857 0.7041858 0.7041859
[211] 0.7041860 0.7041860 0.7041861 0.7041862 0.7041863 0.7041863 0.7041864
[218] 0.7041864 0.7041865 0.7041866 0.7041866 0.7041867 0.7041867 0.7041867
[225] 0.7041868 0.7041868 0.7041869 0.7041869 0.7041869 0.7041870 0.7041870
[232] 0.7041870 0.7041871 0.7041871 0.7041871 0.7041871 0.7041872 0.7041872
[239] 0.7041872 0.7041872 0.7041872 0.7041873 0.7041873 0.7041873 0.7041873
[246] 0.7041873 0.7041873 0.7041874 0.7041874 0.7041874 0.7041874 0.7041874
[253] 0.7041874 0.7041874 0.7041874 0.7041875 0.7041875 0.7041875 0.7041875
[260] 0.7041875 0.7041875 0.7041875 0.7041875 0.7041875 0.7041875 0.7041875
[267] 0.7041875 0.7041875 0.7041875 0.7041875 0.7041875 0.7041876 0.7041876
[274] 0.7041876 0.7041876 0.7041876 0.7041876 0.7041876 0.7041876 0.7041876
[281] 0.7041876 0.7041876 0.7041876 0.7041876 0.7041876 0.7041876 0.7041876
[288] 0.7041876 0.7041876 0.7041876 0.7041876 0.7041876 0.7041876 0.7041876
[295] 0.7041876 0.7041876 0.7041876 0.7041876 0.7041876 0.7041876 0.7041876
For u1
u1
[1] 105.37122 104.28497 103.19437 102.19782 101.34808 100.69007 100.21667 99.88967 99.66755
[10] 99.51631 99.41138 99.33623 99.28022 99.23661 99.20121 99.17139 99.14549 99.12247
[19] 99.10163 99.08255 99.06490 99.04849 99.03315 99.01877 99.00527 98.99257 98.98060
[28] 98.96932 98.95867 98.94862 98.93912 98.93015 98.92166 98.91364 98.90605 98.89886
[37] 98.89206 98.88562 98.87952 98.87374 98.86826 98.86307 98.85815 98.85348 98.84905
[46] 98.84485 98.84087 98.83709 98.83350 98.83009 98.82686 98.82379 98.82087 98.81810
[55] 98.81547 98.81297 98.81059 98.80833 98.80619 98.80415 98.80221 98.80037 98.79861
[64] 98.79695 98.79536 98.79386 98.79243 98.79106 98.78977 98.78853 98.78736 98.78624
[73] 98.78518 98.78417 98.78321 98.78230 98.78143 98.78060 98.77981 98.77906 98.77835
[82] 98.77767 98.77702 98.77640 98.77582 98.77526 98.77473 98.77422 98.77374 98.77329
[91] 98.77285 98.77243 98.77204 98.77166 98.77130 98.77096 98.77064 98.77033 98.77003
[100] 98.76975 98.76949 98.76923 98.76899 98.76876 98.76854 98.76833 98.76813 98.76794
[109] 98.76776 98.76759 98.76743 98.76727 98.76712 98.76698 98.76685 98.76672 98.76660
[118] 98.76648 98.76637 98.76626 98.76616 98.76607 98.76598 98.76589 98.76581 98.76573
[127] 98.76565 98.76558 98.76551 98.76545 98.76539 98.76533 98.76527 98.76522 98.76517
[136] 98.76512 98.76507 98.76503 98.76499 98.76495 98.76491 98.76487 98.76484 98.76481
[145] 98.76478 98.76475 98.76472 98.76469 98.76467 98.76464 98.76462 98.76460 98.76457
[154] 98.76455 98.76454 98.76452 98.76450 98.76448 98.76447 98.76445 98.76444 98.76442
[163] 98.76441 98.76440 98.76439 98.76438 98.76437 98.76436 98.76435 98.76434 98.76433
[172] 98.76432 98.76431 98.76430 98.76430 98.76429 98.76428 98.76428 98.76427 98.76427
[181] 98.76426 98.76425 98.76425 98.76425 98.76424 98.76424 98.76423 98.76423 98.76423
[190] 98.76422 98.76422 98.76422 98.76421 98.76421 98.76421 98.76420 98.76420 98.76420
[199] 98.76420 98.76419 98.76419 98.76419 98.76419 98.76419 98.76419 98.76418 98.76418
[208] 98.76418 98.76418 98.76418 98.76418 98.76418 98.76417 98.76417 98.76417 98.76417
[217] 98.76417 98.76417 98.76417 98.76417 98.76417 98.76417 98.76417 98.76417 98.76416
[226] 98.76416 98.76416 98.76416 98.76416 98.76416 98.76416 98.76416 98.76416 98.76416
[235] 98.76416 98.76416 98.76416 98.76416 98.76416 98.76416 98.76416 98.76416 98.76416
[244] 98.76416 98.76416 98.76416 98.76416 98.76416 98.76416 98.76416 98.76416 98.76416
[253] 98.76415 98.76415 98.76415 98.76415 98.76415 98.76415 98.76415 98.76415 98.76415
[262] 98.76415 98.76415 98.76415 98.76415 98.76415 98.76415 98.76415 98.76415 98.76415
[271] 98.76415 98.76415 98.76415 98.76415 98.76415 98.76415 98.76415 98.76415 98.76415
[280] 98.76415 98.76415 98.76415 98.76415 98.76415 98.76415 98.76415 98.76415 98.76415
[289] 98.76415 98.76415 98.76415 98.76415 98.76415 98.76415 98.76415 98.76415 98.76415
[298] 98.76415 98.76415 98.76415 98.76415
for u2
[1] 133.5728 133.8861 134.2721 134.6356 134.9464 135.1737 135.3096 135.3706 135.3799 135.3572
[11] 135.3157 135.2642 135.2078 135.1499 135.0923 135.0362 134.9820 134.9302 134.8809 134.8340
[21] 134.7895 134.7473 134.7073 134.6695 134.6337 134.5997 134.5676 134.5371 134.5083 134.4809
[31] 134.4550 134.4305 134.4072 134.3851 134.3641 134.3442 134.3254 134.3075 134.2905 134.2743
[41] 134.2590 134.2445 134.2307 134.2175 134.2051 134.1932 134.1820 134.1713 134.1611 134.1515
[51] 134.1423 134.1335 134.1252 134.1173 134.1098 134.1027 134.0959 134.0895 134.0833 134.0775
[61] 134.0719 134.0666 134.0616 134.0568 134.0523 134.0479 134.0438 134.0399 134.0361 134.0326
[71] 134.0292 134.0260 134.0229 134.0200 134.0172 134.0146 134.0121 134.0097 134.0074 134.0052
[81] 134.0032 134.0012 133.9993 133.9975 133.9958 133.9942 133.9927 133.9912 133.9898 133.9885
[91] 133.9872 133.9860 133.9849 133.9838 133.9827 133.9818 133.9808 133.9799 133.9791 133.9782
[101] 133.9775 133.9767 133.9760 133.9754 133.9747 133.9741 133.9735 133.9730 133.9725 133.9720
[111] 133.9715 133.9710 133.9706 133.9702 133.9698 133.9694 133.9691 133.9687 133.9684 133.9681
[121] 133.9678 133.9675 133.9673 133.9670 133.9668 133.9665 133.9663 133.9661 133.9659 133.9657
[131] 133.9655 133.9654 133.9652 133.9651 133.9649 133.9648 133.9646 133.9645 133.9644 133.9643
[141] 133.9642 133.9640 133.9639 133.9639 133.9638 133.9637 133.9636 133.9635 133.9634 133.9634
[151] 133.9633 133.9632 133.9632 133.9631 133.9631 133.9630 133.9630 133.9629 133.9629 133.9628
[161] 133.9628 133.9627 133.9627 133.9627 133.9626 133.9626 133.9626 133.9625 133.9625 133.9625
[171] 133.9625 133.9624 133.9624 133.9624 133.9624 133.9623 133.9623 133.9623 133.9623 133.9623
[181] 133.9623 133.9622 133.9622 133.9622 133.9622 133.9622 133.9622 133.9622 133.9622 133.9621
[191] 133.9621 133.9621 133.9621 133.9621 133.9621 133.9621 133.9621 133.9621 133.9621 133.9621
[201] 133.9621 133.9621 133.9620 133.9620 133.9620 133.9620 133.9620 133.9620 133.9620 133.9620
[211] 133.9620 133.9620 133.9620 133.9620 133.9620 133.9620 133.9620 133.9620 133.9620 133.9620
[221] 133.9620 133.9620 133.9620 133.9620 133.9620 133.9620 133.9620 133.9620 133.9620 133.9620
[231] 133.9620 133.9620 133.9620 133.9620 133.9620 133.9620 133.9620 133.9620 133.9620 133.9620
[241] 133.9620 133.9620 133.9620 133.9620 133.9620 133.9620 133.9620 133.9620 133.9620 133.9620
[251] 133.9620 133.9620 133.9620 133.9620 133.9619 133.9619 133.9619 133.9619 133.9619 133.9619
[261] 133.9619 133.9619 133.9619 133.9619 133.9619 133.9619 133.9619 133.9619 133.9619 133.9619
[271] 133.9619 133.9619 133.9619 133.9619 133.9619 133.9619 133.9619 133.9619 133.9619 133.9619
[281] 133.9619 133.9619 133.9619 133.9619 133.9619 133.9619 133.9619 133.9619 133.9619 133.9619
[291] 133.9619 133.9619 133.9619 133.9619 133.9619 133.9619 133.9619 133.9619 133.9619 133.9619
[301] 133.9619
For sigma1 square:
[1] 220.12489 205.10096 183.53226 160.65836 141.18629 127.45956 118.93847 113.86190 110.77411
[10] 108.80379 107.47407 106.52342 105.80421 105.23069 104.75194 104.33723 103.96770 103.63162
[19] 103.32155 103.03264 102.76167 102.50639 102.26517 102.03679 101.82027 101.61481 101.41971
[28] 101.23438 101.05826 100.89085 100.73171 100.58039 100.43649 100.29964 100.16948 100.04567
[37] 99.92790 99.81586 99.70927 99.60785 99.51136 99.41954 99.33217 99.24902 99.16990
[46] 99.09459 99.02292 98.95471 98.88979 98.82799 98.76917 98.71318 98.65988 98.60914
[55] 98.56084 98.51486 98.47108 98.42941 98.38973 98.35196 98.31599 98.28174 98.24914
[64] 98.21810 98.18854 98.16039 98.13359 98.10807 98.08377 98.06063 98.03860 98.01762
[73] 97.99764 97.97861 97.96049 97.94324 97.92681 97.91116 97.89626 97.88207 97.86855
[82] 97.85568 97.84342 97.83175 97.82063 97.81005 97.79996 97.79036 97.78122 97.77251
[91] 97.76421 97.75631 97.74879 97.74162 97.73479 97.72829 97.72210 97.71621 97.71059
[100] 97.70524 97.70015 97.69530 97.69067 97.68627 97.68208 97.67809 97.67429 97.67067
[109] 97.66722 97.66393 97.66080 97.65782 97.65498 97.65228 97.64970 97.64725 97.64491
[118] 97.64269 97.64057 97.63855 97.63663 97.63480 97.63305 97.63139 97.62981 97.62830
[127] 97.62687 97.62550 97.62420 97.62296 97.62178 97.62065 97.61958 97.61856 97.61759
[136] 97.61666 97.61578 97.61494 97.61414 97.61338 97.61265 97.61196 97.61130 97.61067
[145] 97.61008 97.60951 97.60897 97.60845 97.60796 97.60749 97.60704 97.60662 97.60621
[154] 97.60583 97.60546 97.60511 97.60478 97.60446 97.60416 97.60387 97.60360 97.60334
[163] 97.60309 97.60285 97.60263 97.60241 97.60221 97.60201 97.60183 97.60165 97.60148
[172] 97.60132 97.60117 97.60102 97.60088 97.60075 97.60062 97.60051 97.60039 97.60028
[181] 97.60018 97.60008 97.59999 97.59990 97.59981 97.59973 97.59965 97.59958 97.59951
[190] 97.59944 97.59938 97.59932 97.59926 97.59921 97.59915 97.59910 97.59906 97.59901
[199] 97.59897 97.59893 97.59889 97.59885 97.59881 97.59878 97.59875 97.59872 97.59869
[208] 97.59866 97.59863 97.59861 97.59859 97.59856 97.59854 97.59852 97.59850 97.59848
[217] 97.59846 97.59845 97.59843 97.59841 97.59840 97.59839 97.59837 97.59836 97.59835
[226] 97.59834 97.59832 97.59831 97.59830 97.59829 97.59829 97.59828 97.59827 97.59826
[235] 97.59825 97.59825 97.59824 97.59823 97.59823 97.59822 97.59822 97.59821 97.59821
[244] 97.59820 97.59820 97.59819 97.59819 97.59818 97.59818 97.59818 97.59817 97.59817
[253] 97.59817 97.59816 97.59816 97.59816 97.59816 97.59815 97.59815 97.59815 97.59815
[262] 97.59814 97.59814 97.59814 97.59814 97.59814 97.59814 97.59813 97.59813 97.59813
[271] 97.59813 97.59813 97.59813 97.59813 97.59813 97.59812 97.59812 97.59812 97.59812
[280] 97.59812 97.59812 97.59812 97.59812 97.59812 97.59812 97.59812 97.59812 97.59812
[289] 97.59811 97.59811 97.59811 97.59811 97.59811 97.59811 97.59811 97.59811 97.59811
[298] 97.59811 97.59811 97.59811 97.59811
For sigma2 square
[1] 559.5878 525.6849 505.7973 488.0545 472.0863 458.8586 449.1209 442.7654 439.1009 437.3237
[11] 436.7687 436.9631 437.5961 438.4720 439.4702 440.5184 441.5732 442.6096 443.6137 444.5780
[21] 445.4991 446.3760 447.2092 447.9998 448.7495 449.4602 450.1336 450.7719 451.3767 451.9501
[31] 452.4935 453.0088 453.4975 453.9609 454.4006 454.8177 455.2135 455.5891 455.9457 456.2843
[41] 456.6057 456.9110 457.2010 457.4764 457.7380 457.9866 458.2229 458.4474 458.6608 458.8636
[51] 459.0564 459.2397 459.4140 459.5797 459.7373 459.8872 460.0297 460.1653 460.2943 460.4170
[61] 460.5337 460.6447 460.7503 460.8509 460.9465 461.0375 461.1241 461.2065 461.2850 461.3596
[71] 461.4307 461.4983 461.5627 461.6240 461.6823 461.7378 461.7907 461.8410 461.8889 461.9345
[81] 461.9779 462.0192 462.0586 462.0960 462.1317 462.1657 462.1980 462.2288 462.2581 462.2861
[91] 462.3126 462.3380 462.3621 462.3850 462.4069 462.4277 462.4475 462.4664 462.4844 462.5015
[101] 462.5178 462.5333 462.5481 462.5622 462.5756 462.5883 462.6005 462.6121 462.6231 462.6336
[111] 462.6436 462.6532 462.6622 462.6709 462.6791 462.6869 462.6944 462.7015 462.7083 462.7147
[121] 462.7209 462.7267 462.7323 462.7376 462.7427 462.7475 462.7521 462.7564 462.7606 462.7646
[131] 462.7683 462.7719 462.7754 462.7786 462.7817 462.7847 462.7875 462.7902 462.7927 462.7952
[141] 462.7975 462.7997 462.8018 462.8038 462.8057 462.8075 462.8093 462.8109 462.8125 462.8140
[151] 462.8154 462.8167 462.8180 462.8193 462.8204 462.8216 462.8226 462.8236 462.8246 462.8255
[161] 462.8264 462.8272 462.8280 462.8288 462.8295 462.8302 462.8308 462.8315 462.8321 462.8326
[171] 462.8332 462.8337 462.8342 462.8346 462.8351 462.8355 462.8359 462.8363 462.8366 462.8370
[181] 462.8373 462.8376 462.8379 462.8382 462.8385 462.8387 462.8390 462.8392 462.8394 462.8397
[191] 462.8399 462.8401 462.8402 462.8404 462.8406 462.8407 462.8409 462.8410 462.8412 462.8413
[201] 462.8414 462.8416 462.8417 462.8418 462.8419 462.8420 462.8421 462.8422 462.8422 462.8423
[211] 462.8424 462.8425 462.8425 462.8426 462.8427 462.8427 462.8428 462.8428 462.8429 462.8429
[221] 462.8430 462.8430 462.8431 462.8431 462.8432 462.8432 462.8432 462.8433 462.8433 462.8433
[231] 462.8434 462.8434 462.8434 462.8434 462.8435 462.8435 462.8435 462.8435 462.8435 462.8436
[241] 462.8436 462.8436 462.8436 462.8436 462.8436 462.8437 462.8437 462.8437 462.8437 462.8437
[251] 462.8437 462.8437 462.8437 462.8437 462.8438 462.8438 462.8438 462.8438 462.8438 462.8438
[261] 462.8438 462.8438 462.8438 462.8438 462.8438 462.8438 462.8438 462.8438 462.8438 462.8438
[271] 462.8439 462.8439 462.8439 462.8439 462.8439 462.8439 462.8439 462.8439 462.8439 462.8439
[281] 462.8439 462.8439 462.8439 462.8439 462.8439 462.8439 462.8439 462.8439 462.8439 462.8439
[291] 462.8439 462.8439 462.8439 462.8439 462.8439 462.8439 462.8439 462.8439 462.8439 462.8439
[301] 462.8439
2015年6月17日星期三
EM algorithm to impute missing value and underlying parameter
Genetic linkage data are from Rao, C. R. (1973), Linear Statistical Inference and Its Applications, 2nd Ed., New York: John Wiley & Sons. p368
x<-c(125,18,20,34)
n<-sum(x)
theta.0<-(x[1]-x[2]-x[3]+x[4])/n
theta.0
current_theta<-c()
new_theta<-c()
z11<-c()
z12<-c()
i<-1
current_theta[i]<-theta.0
new_theta[i]<-(x[1]*current_theta[i]+2*x[4]+x[4]*current_theta[i])/(n*current_theta[i]+2*(x[2]+x[3]+x[4]))
new_theta[i]
difference<-abs(new_theta[i]-current_theta[i])
difference
while(difference>0.00000001){
current_theta[i]<-new_theta[i]
new_theta[i+1]<-(x[1]*current_theta[i]+2*x[4]+x[4]*current_theta[i])/(n*current_theta[i]+2*(x[2]+x[3]+x[4]))
difference<-abs(new_theta[i+1]-current_theta[i])
z11[i]<-(0.5/(0.5+0.25*current_theta[i]))*x[1]
z12[i]<-((0.25*current_theta[i])/(0.5+0.25*current_theta[i]))*x[1]
i<-i+1
}
cbind(theta=current_theta,Z11=z11,Z12=z12)
Results: Z11 and Z12 are imputed missing values and thetas are mle estimate of parameter, 7 iterations will get 10^-8 accuracy.
theta Z11 Z12
[1,] 0.6251317 95.23332 29.76668
[2,] 0.6265968 95.18019 29.81981
[3,] 0.6267917 95.17314 29.82686
[4,] 0.6268175 95.17220 29.82780
[5,] 0.6268210 95.17207 29.82793
[6,] 0.6268214 95.17206 29.82794
[7,] 0.6268215 95.17206 29.82794
x<-c(125,18,20,34)
n<-sum(x)
theta.0<-(x[1]-x[2]-x[3]+x[4])/n
theta.0
current_theta<-c()
new_theta<-c()
z11<-c()
z12<-c()
i<-1
current_theta[i]<-theta.0
new_theta[i]<-(x[1]*current_theta[i]+2*x[4]+x[4]*current_theta[i])/(n*current_theta[i]+2*(x[2]+x[3]+x[4]))
new_theta[i]
difference<-abs(new_theta[i]-current_theta[i])
difference
while(difference>0.00000001){
current_theta[i]<-new_theta[i]
new_theta[i+1]<-(x[1]*current_theta[i]+2*x[4]+x[4]*current_theta[i])/(n*current_theta[i]+2*(x[2]+x[3]+x[4]))
difference<-abs(new_theta[i+1]-current_theta[i])
z11[i]<-(0.5/(0.5+0.25*current_theta[i]))*x[1]
z12[i]<-((0.25*current_theta[i])/(0.5+0.25*current_theta[i]))*x[1]
i<-i+1
}
cbind(theta=current_theta,Z11=z11,Z12=z12)
Results: Z11 and Z12 are imputed missing values and thetas are mle estimate of parameter, 7 iterations will get 10^-8 accuracy.
theta Z11 Z12
[1,] 0.6251317 95.23332 29.76668
[2,] 0.6265968 95.18019 29.81981
[3,] 0.6267917 95.17314 29.82686
[4,] 0.6268175 95.17220 29.82780
[5,] 0.6268210 95.17207 29.82793
[6,] 0.6268214 95.17206 29.82794
[7,] 0.6268215 95.17206 29.82794
2015年6月15日星期一
Power and sample size for Binomial distribution
power_n<-function(size){
p0<-51/235 #for calculate the critical region
c<-qbinom(p0,size, 1/6)
p<-seq(0.01,0.4,0.01)
pow<-1-dbinom(c,size,p)+dbinom(size-c,size,p) #the second part is very small can be omitted.
return(pow)
}
p<-seq(0.01,0.4,0.01)
plot(p,power_n(235),type="l",col="red")
lines(p,power_n(300),lty=2,col="blue")
lines(p,power_n(400),lty=3, col="black")
text(0.23,0.93,"Red, sample size =235")
text(0.23,0.94,"Blue, sample size =300")
text(0.24,0.95,"Black, sample size =400")
p0<-51/235 #for calculate the critical region
c<-qbinom(p0,size, 1/6)
p<-seq(0.01,0.4,0.01)
pow<-1-dbinom(c,size,p)+dbinom(size-c,size,p) #the second part is very small can be omitted.
return(pow)
}
p<-seq(0.01,0.4,0.01)
plot(p,power_n(235),type="l",col="red")
lines(p,power_n(300),lty=2,col="blue")
lines(p,power_n(400),lty=3, col="black")
text(0.23,0.93,"Red, sample size =235")
text(0.23,0.94,"Blue, sample size =300")
text(0.24,0.95,"Black, sample size =400")
订阅：
博文 (Atom)