查看: 3637|回复: 3|关注: 0

[已答复] fmincon出现undefined function错误

[复制链接]

新手

5 麦片

财富积分


050


1

主题

2

帖子

0

最佳答案
发表于 2015-12-14 17:21:22 | 显示全部楼层 |阅读模式
第一次发帖希望能得到解决。。。已经调试了两三天了。。。
问题出现在主程序最后一句,样本数据小的时候能顺利的出来结果。但是当我加大样本数量后就会出现
Error using barrier
Objective function is undefined at initial point. Fmincon cannot continue.

Error in fmincon (line 796)
    [X,FVAL,EXITFLAG,OUTPUT,LAMBDA,GRAD,HESSIAN] =
    barrier(funfcn,X,A,B,Aeq,Beq,l,u,confcn,options.HessFcn, ...

Error in fnl_mixed_regress_band_cv (line 77)
    hl_cv = fmincon('fnl_mixed_regress_band_cv_obj',hl_start,[],[],[],[],lb,ub,[],
    options, Y,X,Uc,Ud,W,cases,methods,frequency); % options

以下是主程序

% Modified from Fnl_coeff_mixed1205_cps9005

clear;
tic
%B=100; %number of subsamples in obtaining cross-validated bandwidth
cases=1; % 1 for Epan kernel, 2 for normal kernel
methods=2; % 2 for local constant estimation, not finished yet
          % 1 for local linear estimation:

n_eval=30; %old: 47; % number of evaluation points for plotting purpose (age: 19:65)

warning off  %MATLAB:singularMatrix
    D0=xlsread( 'labor_data.xlsx'); %ID sex edu workexp totincome
    D0=[D0(:,2) D0(:,3) D0(:,4) D0(:,5)]; %sex educ exper income

    [n,p0]=size(D0);   
    Y=log(D0(:,4));
    X1=D0(:,2); %education years
    Z1=zeros(n,1);%no IV
    Uc=D0(:,3);
    Ud=D0(:,1);      
    X=[ones(n,1) X1]; % n x d
    Z=[ones(n,1) Z1]; % n x q
    U=[Uc Ud];


           figure(1);
      %     boxplot(Y,Uc);
      %     xlabel('Experience');
      %     ylabel('Logarithm of wage');
      %     title('Experience-Wage Profile');

      %     figure(2);
      %     boxplot(Y,X1);
      %     xlabel('Years of Schooling');
      %     ylabel('Logarithm of wage');
      %     title('Schooling-Wage Profile');

      %     figure(3);
           subplot 121;
           h=boxplot(Y,Uc);
           xlabel('Experience');
           ylabel('Logarithm of wage');
           title('Experience-Wage Profile');
           subplot 122;
           hl=boxplot(Y,X1);
           xlabel('Years of Schooling');
           ylabel('Logarithm of wage');
           title('Schooling-Wage Profile');
          % figure(10);
          % plot(Uc,Y);
          % xlabel('Experience');       %  add axis labels and plot title
          % ylabel('logarithm of wage');
          % title('Experience-Wage Profile');


    [n,d]=size(X);
    [n,p]=size(Uc);
    [n,q]=size(Ud);
    [n,qz]=size(Z);

   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Functional coefficient estimator      

    W=ones(n,p);  

       for ii=1:p
           aa=Uc(:,ii);
           Ucmean=mean(aa);
           Ucstd=std(aa);
           W(:,ii)= 1*( abs(aa-Ucmean)<=2*Ucstd );
       end
       W=prod(W,2); % n x 1      

       xu=[X(:,2:d) Uc];
       for ii=1:(d-1+p)
           aa=xu(:,ii);
           Ucmean=mean(aa);
           Ucstd=std(aa);
       end      

    % Our nonparametric method
    % rule of thumb
  %  h_msu_tilde=2.37*std(Uc)*n^(-1/(p+4));  
  %  lambda_msu_tilde=1*std(Uc)*n^(-2/(p+4))*ones(1,q);   
  %  h_msu_hat=2.37*std(Uc)*n^(-1/(p+4));  
  %  lambda_msu_hat=1*std(Uc)*n^(-2/(p+4))*ones(1,q);   

    % For Su, Chen and Ullah's method
    [h_scu,lambda_scu] = fnl_mixed_regress_band_cv(Y,X,Uc,Ud,W,cases,methods,0);% for CV (h,lambda)


以下是fnl_mixed_regress_band_cv
function [hcv, lambdacv]=fnl_mixed_regress_band_cv(Y,X,Uc,Ud,W,cases,methods,frequency)
% Purpose:  Least squares cross validation for nonparametric regression with mixed regressors
%Inputs:
% Case (1) For mixed regressors in functional coefficients:
% Y: dependent variable  (n x 1)
% X: regressor ( n x d)
% Uc: continuous independent variables  (n x p)
% Ud: discrete independent variables  (n x q)
% W: Weight function ( n x 1)
% cases=1: 2nd order Epan kernel,
% cases=2: 2nd Gaussian kernel,
% cases=3: 4th order Epan kernel,
% cases=4: 4th order Gaussian kernel
% methods= 1: local constant estimation % to be done
% methods= 2: local linear estimation
% frequency=0: cross-validate both h and lambda, hl: (p+q)-vector
% frequency=1: cross-validate only h, set lambda=0, hl: p-vector


% Case (2) For pure continuous regressors in functional coefficients:
% Ud=[]; %empty matrix

% Output:
% hcv: cross validated bandwidth (1xp) for continuous variables
% lambdacv: cross validated bandwidth (1xq) for discrete variables
% Code written by Liangjun Su
% Guanghua School of Management, Peking University
% Beijing 100871, China
% Email: lsu%gsm.pku.edu.cn
% Last modified: 3/8/06

if isempty(Ud)==0
    [n,d]=size(X);
    [n,p]=size(Uc); [n,q]=size(Ud);
    %n: #of rows=number of observations

    if cases>4; error('cases should be less than 4'); end
    if methods>2 error('method should be 1 or 2: 1 for local constant estiamtion and 2 for local linear estimation');end

    if n<=1000;
        options=optimset('maxiter',1000,'tolx',2e-5,'tolfun', 2e-5,...
             'display','off','LargeScale'  , 'off');  %the first three numbers control the speed for each repetition
    else;
        options=optimset('maxiter',500,'tolx', 2e-5,'tolfun', 2e-5,...
             'display','off','LargeScale'  , 'off');  %the first three numbers control the speed for each repetition
    end        

    % Bound for h
     if cases==1;     lbc=0.5*ones(1,p); ubc=5*ones(1,p);
     elseif cases==2; lbc=0.1*ones(1,p); ubc=5*ones(1,p);
     end
    % Bound for lambda
    lbd=0*zeros(1,q);
    ubd=ones(1,q);

    % for s=1:q;
    %     temp=unique(Ud(:,s)); % count number of categories in Ud
    %     cs=length(temp); %L0(1,s)=length(temp);
    %     ubd(1,s)=(cs-1)/cs; % upper bound for lambda_s
    % end


    % assign starting values
    h_start=ones(1,p);
    if frequency==0
          lambda_start=0.5*n^(-2/(4+p))*ones(1,q);
          hl_start=[h_start lambda_start];
          lb=[lbc lbd]; ub=[ubc ubd];
    elseif frequency==1
          % lambda_start=0.5*n^(-2/(4+p))*ones(1,q);
          hl_start=h_start;
          lb=lbc; ub=ubc;
    else
    end


    hl_cv = fmincon('fnl_mixed_regress_band_cv_obj',hl_start,[],[],[],[],lb,ub,[], options, Y,X,Uc,Ud,W,cases,methods,frequency); % options
    %size(hl_cv)
    if cases==1|cases==2
        hcv=hl_cv(1:p).*std(Uc)*n^(-1/(4+p)); % for 2nd order kernel: 1 x p
    elseif cases==3|cases==4
        hcv=hl_cv(1:p).*std(Uc)*n^(-1/(8+p)); % for 4th order kernel: 1 x p
    else
    end
    if frequency==0;
         lambdacv=hl_cv(p+1:p+q);
          %lambdacv=lambdacv*0.1*n^(-2/(p+4));
    else
         lambdacv=zeros(1,q);
    end

elseif isempty(Ud)==1

    [n,d]=size(X);
    [n,p]=size(Uc);
    %n: #of rows=number of observations

    if cases>4; error('cases should be less than 4'); end
    if methods>2 error('method should be 1 or 2: 1 for local constant estiamtion and 2 for local linear estimation');end

    if n<=1000;
        options=optimset('maxiter',1000,'tolx',2e-5,'tolfun', 2e-5,...
             'display','off','LargeScale'  , 'off');  %the first three numbers control the speed for each repetition
    else;
        options=optimset('maxiter',500,'tolx', 2e-5,'tolfun', 2e-5,...
             'display','off','LargeScale'  , 'off');  %the first three numbers control the speed for each repetition
    end        

    % Bound for h
     if cases==1;     lbc=0.5*ones(1,p); ubc=5*ones(1,p);
     elseif cases==2; lbc=0.1*ones(1,p); ubc=5*ones(1,p);
     end
%     % Bound for lambda
%     lbd=0*zeros(1,q);
%     ubd=ones(1,q);
%
%     lb=[lbc lbd]; ub=[ubc ubd];

    % assign starting values
    h_start=ones(1,p);


    h1_cv = fminunc('fnl_mixed_regress_band_cv_obj', h_start, options, Y,X,Uc,Ud,W,cases,methods); % options

    if cases==1|cases==2
        hcv=h1_cv(1:p).*std(Uc)*n^(-1/(4+p)); % for 2nd order kernel: 1 x p
    elseif cases==3|cases==4
        hcv=hl_cv(1:p).*std(Uc)*n^(-1/(8+p)); % for 2nd order kernel: 1 x p
    else
    end
    lambdacv=[];

else

end



以下是fnl_mixed_regress_band_cv_obj
function [hcv, lambdacv]=fnl_mixed_regress_band_cv(Y,X,Uc,Ud,W,cases,methods,frequency)
% Purpose:  Least squares cross validation for nonparametric regression with mixed regressors
%Inputs:
% Case (1) For mixed regressors in functional coefficients:
% Y: dependent variable  (n x 1)
% X: regressor ( n x d)
% Uc: continuous independent variables  (n x p)
% Ud: discrete independent variables  (n x q)
% W: Weight function ( n x 1)
% cases=1: 2nd order Epan kernel,
% cases=2: 2nd Gaussian kernel,
% cases=3: 4th order Epan kernel,
% cases=4: 4th order Gaussian kernel
% methods= 1: local constant estimation % to be done
% methods= 2: local linear estimation
% frequency=0: cross-validate both h and lambda, hl: (p+q)-vector
% frequency=1: cross-validate only h, set lambda=0, hl: p-vector


% Case (2) For pure continuous regressors in functional coefficients:
% Ud=[]; %empty matrix

% Output:
% hcv: cross validated bandwidth (1xp) for continuous variables
% lambdacv: cross validated bandwidth (1xq) for discrete variables
% Code written by Liangjun Su
% Guanghua School of Management, Peking University
% Beijing 100871, China
% Email: lsu%gsm.pku.edu.cn
% Last modified: 3/8/06

if isempty(Ud)==0
    [n,d]=size(X);
    [n,p]=size(Uc); [n,q]=size(Ud);
    %n: #of rows=number of observations

    if cases>4; error('cases should be less than 4'); end
    if methods>2 error('method should be 1 or 2: 1 for local constant estiamtion and 2 for local linear estimation');end

    if n<=1000;
        options=optimset('maxiter',1000,'tolx',2e-5,'tolfun', 2e-5,...
             'display','off','LargeScale'  , 'off');  %the first three numbers control the speed for each repetition
    else;
        options=optimset('maxiter',500,'tolx', 2e-5,'tolfun', 2e-5,...
             'display','off','LargeScale'  , 'off');  %the first three numbers control the speed for each repetition
    end        

    % Bound for h
     if cases==1;     lbc=0.5*ones(1,p); ubc=5*ones(1,p);
     elseif cases==2; lbc=0.1*ones(1,p); ubc=5*ones(1,p);
     end
    % Bound for lambda
    lbd=0*zeros(1,q);
    ubd=ones(1,q);

    % for s=1:q;
    %     temp=unique(Ud(:,s)); % count number of categories in Ud
    %     cs=length(temp); %L0(1,s)=length(temp);
    %     ubd(1,s)=(cs-1)/cs; % upper bound for lambda_s
    % end


    % assign starting values
    h_start=ones(1,p);
    if frequency==0
          lambda_start=0.5*n^(-2/(4+p))*ones(1,q);
          hl_start=[h_start lambda_start];
          lb=[lbc lbd]; ub=[ubc ubd];
    elseif frequency==1
          % lambda_start=0.5*n^(-2/(4+p))*ones(1,q);
          hl_start=h_start;
          lb=lbc; ub=ubc;
    else
    end


    hl_cv = fmincon('fnl_mixed_regress_band_cv_obj',hl_start,[],[],[],[],lb,ub,[], options, Y,X,Uc,Ud,W,cases,methods,frequency); % options
    %size(hl_cv)
    if cases==1|cases==2
        hcv=hl_cv(1:p).*std(Uc)*n^(-1/(4+p)); % for 2nd order kernel: 1 x p
    elseif cases==3|cases==4
        hcv=hl_cv(1:p).*std(Uc)*n^(-1/(8+p)); % for 4th order kernel: 1 x p
    else
    end
    if frequency==0;
         lambdacv=hl_cv(p+1:p+q);
          %lambdacv=lambdacv*0.1*n^(-2/(p+4));
    else
         lambdacv=zeros(1,q);
    end

elseif isempty(Ud)==1

    [n,d]=size(X);
    [n,p]=size(Uc);
    %n: #of rows=number of observations

    if cases>4; error('cases should be less than 4'); end
    if methods>2 error('method should be 1 or 2: 1 for local constant estiamtion and 2 for local linear estimation');end

    if n<=1000;
        options=optimset('maxiter',1000,'tolx',2e-5,'tolfun', 2e-5,...
             'display','off','LargeScale'  , 'off');  %the first three numbers control the speed for each repetition
    else;
        options=optimset('maxiter',500,'tolx', 2e-5,'tolfun', 2e-5,...
             'display','off','LargeScale'  , 'off');  %the first three numbers control the speed for each repetition
    end        

    % Bound for h
     if cases==1;     lbc=0.5*ones(1,p); ubc=5*ones(1,p);
     elseif cases==2; lbc=0.1*ones(1,p); ubc=5*ones(1,p);
     end
%     % Bound for lambda
%     lbd=0*zeros(1,q);
%     ubd=ones(1,q);
%
%     lb=[lbc lbd]; ub=[ubc ubd];

    % assign starting values
    h_start=ones(1,p);


    h1_cv = fminunc('fnl_mixed_regress_band_cv_obj', h_start, options, Y,X,Uc,Ud,W,cases,methods); % options

    if cases==1|cases==2
        hcv=h1_cv(1:p).*std(Uc)*n^(-1/(4+p)); % for 2nd order kernel: 1 x p
    elseif cases==3|cases==4
        hcv=hl_cv(1:p).*std(Uc)*n^(-1/(8+p)); % for 2nd order kernel: 1 x p
    else
    end
    lambdacv=[];

else
end


以下是数据
sex  edu  workexp   income
2        12        20        9600
1        9        41        8400
1        15        46        11851
2        4        30        4761
1        11        31        24684
2        11        32        20067
2        13        2        7015
1        9        17        5811
2        12        21        6206
1        12        36        18057
2        9        28        9559
2        15        1        0
1        12        32        10255
2        12        31        10141
1        13        2        5056
1        11        35        7732
2        9        32        7057
1        12        1        5014
1        11        0        91
1        13        21        11248
2        13        21        10243
2        5        0        8
1        13        31        7040
2        9        26        6242
1        6        29        6282
2        0        15        3223
1        9        31        19122
2        5        23        5040
1        9        0        5360
1        17        35        21027
2        15        30        11835
2        6        30        11653
2        14        2        5414
1        14        15        6150
2        14        13        15950
1        10        11        9073
2        10        11        9430
1        10        16        5828
2        10        18        5532
1        14        28        20860
2        15        24        9183
2        7        30        9350
1        8        42        3530
2        9        37        5732
1        12        37        7416
1        5        36        4720
2        5        28        4720
1        4        45        6759
2        0        0        75
2        11        28        7702
1        4        0        575
1        8        44        11785
2        8        37        3800
1        11        9        9283
1        6        46        4721
2        8        35        4353
1        13        8        6291
1        9        27        6551
2        9        23        7902
1        16        38        13350
2        11        33        5930
1        16        10        7408
2        16        10        5097
1        9        15        7327
2        12        10        8039
2        0        0        8
1        15        41        9767
2        11        35        8849
1        15        0        30
1        11        27        5868
2        8        26        5522
2        11        25        5962
1        11        26        6370
1        11        26        6533
2        11        23        6352
1        12        25        8181
2        14        27        8374
1        17        29        9472
2        15        28        10102
2        12        32        5030
1        14        38        4795
2        6        37        4965
1        5        41        7780
2        5        37        4998
1        8        45        5383
1        11        28        9020
2        14        25        11730
1        11        9        9098
2        12        6        11700
2        0        0        0
1        12        38        8591
2        9        30        7406
2        11        37        5022
1        7        39        8173
1        11        2        3078
2        4        30        3990
1        15        2        7770
2        9        32        9034
1        8        39        4135
1        8        40        9388
1        11        15        8966
2        11        14        7820
1        12        6        16020
2        11        7        9181
1        4        36        6299
2        5        29        4830
2        9        25        10008
1        14        24        8224
1        7        0        463
1        9        31        4809
2        2        10        2084
1        9        22        3718
2        12        15        3890
1        3        0        90
1        9        31        16503
2        12        23        3350
2        10        0        90
1        9        33        7941
2        12        5        4287
2        9        27        1978
1        12        15        5949
2        12        15        7994
2        2        8        5373
1        8        34        7143
2        9        28        5081
1        10        35        4783
1        12        3        1505
1        12        26        13189
2        9        26        10393
2        8        0        36
1        9        27        9592
2        15        28        9357
2        10        0        90
2        14        28        11933
1        12        27        9046
1        0        0        12711
2        0        0        7010
1        16        26        19927
2        16        26        7380
1        14        0        83
2        14        0        545
2        12        37        15790
1        10        34        15037
1        12        0        23962
2        12        23        8981
1        12        26        8720
2        15        10        8914
1        15        13        6845



帖子有点长。。。谢谢大家阅读。。。。。

新手

5 麦片

财富积分


050


1

主题

2

帖子

0

最佳答案
 楼主| 发表于 2015-12-14 17:22:28 | 显示全部楼层
如果取前10个样本数据就不会出现错误,但是我用50个样本数据就会出现错误。

新手

5 麦片

财富积分


050


0

主题

2

帖子

0

最佳答案
发表于 2018-1-25 10:51:49 | 显示全部楼层
flora_zhang 发表于 2015-12-14 17:22
如果取前10个样本数据就不会出现错误,但是我用50个样本数据就会出现错误。 ...

楼主您好!请问问题解决了吗?我也遇到了同样的问题,希望能交流一下!谢谢~

新手

10 麦片

财富积分


050


4

主题

21

帖子

0

最佳答案
发表于 2019-11-7 10:29:45 | 显示全部楼层
我用优化求解也是这个问题,有大佬解决了吗?
您需要登录后才可以回帖 登录 | 注册

本版积分规则

关闭

站长推荐上一条 /3 下一条

快速回复 返回顶部 返回列表