各位大神好，小弟想估计如下的似然函数，但是遇到题中所描述的问题，想请问下怎么解决，看起来应该是说矩阵的维度不对，我仔细研究了一下，因为是有两组待估参数beta，我不清楚这个在R里面怎么去表示，敬请指教
代码中有其它问题，也烦请指出，谢谢！

f(yit∣xit)=pσ1ϕ(yit−xitβ1σ1)+1−pσ2ϕ(yit−xitβ2σ2)f(y_{it}|x_{it})=\frac{p}{\sigma _{1}}\phi \left ( \frac{y_{it}- x_{it}\beta _{1}}{\sigma _{1}}\right )+\frac{1-p}{\sigma _{2}}\phi \left ( \frac{y_{it}- x_{it}\beta _{2}}{\sigma _{2}}\right )

代码如下：

slope1 <- -.3;slope2 <- .3;slope3 <- 1.8; slope4 <- 0.5;intercept1 <- 1.5    
age <- sample(seq(-2,2,len=201), 400)
grade <- sample(seq(-2,2,len=201), 400)
not_smsa <- sample(seq(-2,2,len=201), 400)
wage <- intercept1 + slope1*age +slope2*grade + slope3*not_smsa + rnorm(length(age),0,.15) 

X <- cbind(1, age , grade , not_smsa )
Y <- wage
N<- length(y)
mydata <- cbind.data.frame(X,Y)

ans <- lm(ln_wage ~ age + grade + not_smsa + south, 
          data = xtlogit_data)
vi <- c(coef(ans) , coef(ans) , 0.5)

fmm.nll <- function (beta) {
  
  mu1 <- X %*% beta[1]
  sigma1 <- sqrt((t((Y-X %*% beta[1]))%*%(Y-X %*% beta[1]))/N)
  
  mu2 <- X %*% beta[2]
  sigma2 <- sqrt((t((Y-X %*% beta[2]))%*%(Y-X %*% beta[2]))/N)

  c1 <- dnorm(Y,mu1,sigma1)
  c2 <- dnorm(Y,mu1,sigma1)
  pi <- exp(mu1)/(exp(mu1)+exp(mu2))
  tot=pi*c1+(1-pi)*(c2)
  -sum(log(tot),na.rm=T)
}


fit <- optim(vi, fmm.nll , method = "BFGS", hessian = TRUE)
fit$par

    zhangice001

    可能是我之前没有解释清楚。你的代码里有些比较明显的问题，主要是用了矩阵乘法；另外 beta1 和 beta2 要分别定义。可能还有其他小问题，我没有细看，你先参考下。

    slope1 <- -.3;slope2 <- .3;slope3 <- 1.8; slope4 <- 0.5;intercept1 <- 1.5    
    
age <- sample(seq(-2,2,len=201), 400) # sample size is larger than the vector length without replacing
grade <- sample(seq(-2,2,len=201), 400) # sample size is larger than the vector length without replacing
not_smsa <- sample(seq(-2,2,len=201), 400) # sample size is larger than the vector length without replacing
wage <- intercept1 + slope1*age +slope2*grade + slope3*not_smsa + rnorm(length(age),0,.15) 

X <- cbind(1, age , grade , not_smsa )
Y <- wage
N<- length(y) # here y is not defined
mydata <- cbind.data.frame(X,Y)

ans <- lm(ln_wage ~ age + grade + not_smsa + south, 
          data = xtlogit_data) # xtlogit_data not found
vi <- c(coef(ans) , coef(ans) , 0.5)

fmm.nll <- function (beta1, beta2) { # You need to define beta 1 and beta2 seperately
  
    
  mu1 <- X * beta1 # In the answer, I said %*% was for matrix multiplication, please use * here
  sigma1 <- sqrt((t((Y-X * beta1)) * (Y-X * beta1))/N) #
  
    
  mu2 <- X * beta2
  sigma2 <- sqrt((t((Y-X * beta2)) * (Y-X * beta2))/N)
  
    
  c1 <- dnorm(Y,mu1,sigma1)
  c2 <- dnorm(Y,mu1,sigma1)
  pi <- exp(mu1)/(exp(mu1)+exp(mu2))
  tot=pi*c1+(1-pi)*(c2)
  -sum(log(tot),na.rm=T)
}


fit <- optim(vi, fmm.nll , method = "BFGS", hessian = TRUE)
fit$par
```

        Liechi 谢谢您的帮助和指导。十分抱歉，早上将代码中数据替换成模拟数据时犯了一些低级错误。我看了下您的修改，重新跑了一下，结果显示“non-conformable arrays In addition: Warning messages:”，似乎还是矩阵维度的问题？但是这里模拟生成的应该是没问题的啊？麻烦您了。

        slope1 <- -.3;slope2 <- .3;slope3 <- 1.8; slope4 <- 0.5;intercept1 <- 1.5    
age <- sample(seq(-2,2,len=201), 40) 
grade <- sample(seq(-2,2,len=201), 40) 
not_smsa <- sample(seq(-2,2,len=201), 40) 
wage <- intercept1 + slope1*age +slope2*grade + slope3*not_smsa + rnorm(length(age),0,.15) 

X <- cbind(1, age , grade , not_smsa )
Y <- wage
N<- length(Y) 
mydata <- cbind.data.frame(X,Y)

ans <- lm(wage ~ age + grade + not_smsa , 
          data = mydata) 
vi <- c(coef(ans) , coef(ans) , 0.5)

fmm.nll <- function (beta1, beta2) { # You need to define beta 1 and beta2 seperately
  
  
  mu1 <- X * beta1 # In the answer, I said %*% was for matrix multiplication, please use * here
  sigma1 <- sqrt((t((Y-X * beta1)) * (Y-X * beta1))/N) #
  
  
  mu2 <- X * beta2
  sigma2 <- sqrt((t((Y-X * beta2)) * (Y-X * beta2))/N)
  
  
  c1 <- dnorm(Y,mu1,sigma1)
  c2 <- dnorm(Y,mu2,sigma2)
  pi <- exp(mu1)/(exp(mu1)+exp(mu2))
  tot=pi*c1+(1-pi)*(c2)
  -sum(log(tot),na.rm=T)
}


fit <- optim(vi, fmm.nll , method = "BFGS", hessian = TRUE)
fit$par

            zhangice001

            sigma1 <- sqrt((t((Y-X * beta1)) * (Y-X * beta1))/N) 这一句有些问题：

            X 是一个矩阵，Y 是个向量，直接加减前确定一下那是你想要的吗？如果你希望求标准差的话，公式用错了，t() 转置了啥？你可以从矩阵里选择你需要的列来加减，例如，你可以用 X[ , "age"] 或者 X[ , 2] 来获得 age 那列的数据。

            你好像对 R 的基本操作不熟悉，先花点时间看看基础或许会为以后省下很多时间。别的我没有细看，请你检查一下还有没有类似的问题。

                Liechi 早上好，我仔细研究了一下，其实解决的方法很简单，function里面只要定义一个beta就好了，然后beta1改成beta[1:4]，然后beta2改成[5:8]，还是用矩阵的形式来算就好了