各位好，我用R手写了一个finite mixture model的极大似然估计，但是算出来结果跟stata里面手写的极大似然估计的结果不对，原因出在在R里面，我不知道如何去写每个class的概率的限制条件，烦请指教。

R里面的代码如下：


library(haven)
library(summarytools)
library(matrixStats)

setwd("C:\\Users\\ThinkPad\\Desktop\\MSL for R")
xtprobit_data <- read_dta("xtprobit_data.dta")

y <- unlist(xtprobit_data["ln_wage"]) 
X <- cbind(1, as.matrix(xtprobit_data[c("age", "grade", "not_smsa", "south")]))
#initial values
anso <- lm(ln_wage ~ age + grade + not_smsa + south, 
           data = xtprobit_data)
vi <- c(coef(anso),0.22,0.1,0.1,0.4,0.3,0.4)
#function
fmm <- function(beta) {

  mu1 <- c(X %*% beta[1:5])
  mu2 <- c(X %*% beta[6:10])
  
  p1 <- beta[11]
  p2 <- 1-p1
  
  llk <- p1*dnorm(y,mu1)+p2*dnorm(y,mu2)
  -sum(log(llk),na.rm=T)
}
fit <- optim(vi,fmm , method = "BFGS", hessian = TRUE)
fit$par

在stata里面，对p1、p2只要加这样一个限制条件即可

`p1'=exp(`lp')/(1+exp(`lp'))
`p2'=1-`p1'

stata里面的正确结果和数据在下面的链接里面：链接：https://pan.baidu.com/s/1PwfO4roaiW1CCUuvLB7xBg
提取码：k2h1
复制这段内容后打开百度网盘手机App，操作更方便哦

请问stata这里的限制条件如何在R当中写出来，谢谢！

有没有大佬来救救孩子啊

    6 天 后

    Ihavenothing 您好！谢谢您的指导，我这么写了之后，结果还是跟stata出来的不一样，第一个component和第二个component的结果极为类似，而且p1的值求出来也是为负的，而实际上p1应该是为正的，我将方法改成了“L-BFGS-B”然后加入取值范围限制之后也不行，请问是不是我程序有其它问题？谢谢

      Ihavenothing

      我使用flexmix估计的结果和这个程序的结果完全不一样，将代码粘贴如下，希望您能够给予指导，非常感谢：

      #prepare data
slope1 <- -.3;slope2 <- .3;slope3 <- 1.8; slope4 <- 0.5;intercept1 <- 1.5    
age <- sample(seq(18,60,len=401), 200) 
grade <- sample(seq(0,100,len=401), 200) 
not_smsa <- sample(seq(-2,2,len=401), 200) 
unemployment <- rnorm(200,mean=0,sd=1)
wage <- intercept1 + slope1*age +slope2*grade + slope3*not_smsa + rnorm(length(age),0,.15) 
y <- wage 
X  <- cbind(1, age , grade , not_smsa)
mydata <- cbind.data.frame(X,y)
anso <- lm(wage ~ age + grade + not_smsa, 
           data = mydata)
vi <- c(coef(anso),0.01,0.02,0.03,0.04,0.1)
#function
fmm <- function(beta) {
  
  mu1 <- c(X %*% beta[1:4])
  mu2 <- c(X %*% beta[5:8])
  
  p1 <- 1 / (1 + exp(-beta[9]))
  p2 <- 1-p1
  
  llk <- p1*dnorm(y,mu1)+p2*dnorm(y,mu2)
  -sum(log(llk),na.rm=T)
}

fit <- optim(vi,fmm , method = "BFGS", control = list(maxit=50000), hessian = TRUE)
fit$par

library(flexmix)
flexfit <- flexmix(wage ~ age + grade + not_smsa, data = mydata, k = 2)
flexfit$par
c1 <- parameters(flexfit,component=1)
c2 <- parameters(flexfit, component=2)

          zhangice001 原因很简单，你在似然函数中假定误差方差是1，而在 flexmix 中方差是未知参数。改成下面这样就行了。

          set.seed(123)
slope = c(-0.3, 0.3, 1.8)
intercept = 1.5
sigma = 0.15

age = sample(seq(18, 60, len = 401), 200) 
grade = sample(seq(0, 100, len = 401), 200) 
not_smsa = sample(seq(-2, 2, len = 401), 200)

X = cbind(1, age, grade, not_smsa)
y = wage = c(X %*% c(intercept, slope)) + rnorm(length(age), 0, sigma)

mydata = cbind.data.frame(X, y)
init = lm(wage ~ age + grade + not_smsa, data = mydata)
par_init = c(coef(init), 0, coef(init), 0, 0)

fmm = function(pars) {
    beta1 = pars[1:4]
    sigma1 = log(1 + exp(pars[5]))   ## always positive
    beta2 = pars[6:9]
    sigma2 = log(1 + exp(pars[10]))  ## always positive
    p1 = 1 / (1 + exp(-pars[11]))    ## always in (0, 1)
    
    mu1 = c(X %*% beta1)
    mu2 = c(X %*% beta2)
    
    llk = p1 * dnorm(y, mu1, sigma1) + (1 - p1) * dnorm(y, mu2, sigma2)
    -sum(log(llk))
}

fit = optim(par_init, fmm , method = "BFGS", hessian = TRUE)
fit$par

          还有就是一个合适的初值非常重要。

              zhangice001 其实我很早就看到这个帖子了，但是发现统计模型既没公式又没可重复的小例子，只有一堆代码，就果断放弃了。如果能按照论坛左上角的新手须知添加必要的公式、数据和代码，肯定早有人回复了！

                  Ihavenothing 好的，谢谢大佬！主要是我个人刚刚从stata编程转过来同R，两者的思维逻辑差异好大啊

                    Cloud2016 抱歉抱歉，之前也在本论坛发过类似的帖子，当时是带上了公式，这次偷懒就没有带上