Dear statalist:
I'm sorry to bother you about my question.I want to perform a finite mixture model which is consisted of two ols and their constraints are same.Although stata can do this easily,my further research needs numerical integration and discretization in matlab.But I have failed to do this in matlab.I wonder if you can give some help in writing this fmm in matlab.My data is attached.
And the code in stata is in the follwing:

----------------------- copy starting from the next line -----------------------
Code:
* Example generated by -dataex-. To install: ssc install dataex
*fmm is estimated by MLE without EM
capture program drop fmmreg
program fmmreg
qui {
   args lnf xb1 lns1 xb2 lns2 lp 
   tempvar f1 f2
   gen double `f1'=normalden($ML_y1,`xb1',exp(`lns1'))
   gen double `f2'=normalden($ML_y1,`xb2',exp(`lns2'))
   tempvar p
   gen double `p'=exp(`lp')/(1+exp(`lp'))
   replace `lnf'=ln(`p'*`f1'+(1-`p')*`f2') 
}
end


constraint 1 [xb1]age=[xb2]age
constraint 2 [xb1]education=[xb2]education
constraint 3 [xb1]married=[xb2]married
constraint 4 [xb1]children=[xb2]children
constraint 5 [xb1]county=[xb2]county

ml model lf fmmreg (xb1:wage=age education married children county) (lns1:) (xb2:wage=age education married children county) (lns2:) (lp:),constraints(1/5)
ml max


end
------------------ copy up to and including the previous line ------------------

my code in matlab is in the following:

----------------------- copy starting from the next line -----------------------
Code:
* Example generated by -dataex-. To install: ssc install dataex
data=[lnwage,age,children,county,education,married];
index=sum(isnan(data),2);
data=data(~index,:);
y = data(:,1); 
X = data(:, 2:6); 
[b,std]=ols(y,X,0);
para0=[b;0;0;std;std;0.5];
[para,std]=my_mle(@(para)gmmreg2(para,0,X,y),para0);


function [b,std] = ols( y,X,constant )

if constant == 1
    X = [ones(length(y),1),X];
end
[N,K]=size(X);
b = inv(X'*X)*X'*y;

yhat = X*b;
residul = y - yhat;
std = sqrt((residul'*residul)/(N-K));    
end


function [ para,std] = my_mle(fun,para0)
% para0,initial values
% fv,maximized likelihood value
options=optimset('MaxFunEvals',5000,'MaxIter',5000);
para0=para0(:);
 [para]=fminsearch(fun,para0,options);

d= my_numerical_derivative(fun,para);
std=sqrt(diag(pinv(d'*d)));
end



function f = gmmreg2(para,num,X,y )
b=para(1:5);
mu1=para(6);
mu2=para(7);
sig1=para(8);
sig2=para(9);
p=para(10);
eta=y-X*b;
pdf=p*normpdf(eta,mu1,sig1)+(1-p)*normpdf(eta,mu2,sig2);
pdf=max(pdf,eps);
if num==1
    f=log(pdf);
    else
    f=-sum(log(pdf));
end


function f = my_numerical_derivative(fun,para)
n = length(para);
for i=1:n
    a=zeros(n,1);
    a(i)=max(para(1)*1e-6,1e-6);
    y1(:,i)=feval(fun,para+a);
    y2(:,i)=feval(fun,para-a);
    f(:,i)=(y1(:,i)-y2(:,i))/2/a(i);
end




end
------------------ copy up to and including the previous line ------------------



Thanks,I am looking forward to your reply