声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 8272|回复: 10

[混合编程] index out of bounds because numel(y)=1.

[复制链接]
发表于 2011-3-20 19:42 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

x
%#####Optimal Experimental Design based on the Fisher Information Matrix####
%#####The NF-kB model with 24 states and 64 parameters######################
clc;
close all;
clear all;
format long;
parameters_hoff; % parametrization of ODE's and simulation times


IKK_normal=0.9;  %the initial concentration of IKK

para=[3 4 5];  %the id of parameters analysed

  for k=1:100
      
        Y01(1)=IKK_normal*k/10;
        
        IKK(k)=Y01(1);
        
        tspan=[t0:tau:t0+tw];
        
        [T,Y]=ode15s(@model_hoff,tspan,Y01,[],par,Source);
      
    %%The local sensitivity analysis method is used to gain the absolute sensitivity matrix
      
        error=10;%input gain value[0-100]
        
        disp('please wait......')
        
        abso_s=[];  %the absolute sensitivity matrix
        
    %%the local sensitivity analysis using the infite difference method%%%%%%%%   
     
     for j=1:length(para)
         
          i=para(j);  
         
          %parameter i is decreased with error%
            par(i,1)=par(i,1)*(1-error/100);
          %############ IKK on ####################################
            [T1,Y1]=ode15s(@model_hoff,tspan,Y01,[],par,Source);
          %parameter i is returned to its original value
            par(i,1)=par(i,1)/(1-error/100);
            
          %parameter i is increased with 1error%
            par(i,1)=par(i,1)*(1+error/100);
         %############ IKK on ####################################
            [T2,Y2]=ode15s(@model_hoff,tspan,Y01,[],par,Source);
         %parameter i is returned to its original value
            par(i,1)=par(i,1)/(1+error/100);   
            
         %the absolute sensitivity S(i,j)=[x(p+error%*p)-x(p-error%*p)]/(2*error%*p)
            S=(Y2-Y1)/(2*par(i,1)*error/100);
                 
            abso_s=[abso_s S(:,9)]; %the absolute sensitivity matrix of the system output states x9 w.r.t the i-th parameter        
     end
      
          relative_error=0.05;  %the relative error of the measure value
          absolute_error=0.001; %the absolute error of the measure value

          cova_v=[]; %the error covariance matrix
          v=[]; %the diagnoal elements of the error covariance matrix
         
          for j=9
             for i=1:max(size(Y))
                 sigma(:,i)=relative_error*Y(i,j)+absolute_error;
             end
             v=[v sigma];
          end
          cova_v=diag(v.^2);
  
          %the Fisher information matrix FIM
            FIM=abso_s'*inv(cova_v)*abso_s;
            
            [vec,lam]=eig(FIM); %to calculate the eigenvalue of FIM
            
            for j=1:length(para)
                lambda(j)=lam(j,j); %the eigenvalue lamda of FIM
            end
            
          %Optimal experimental design w.r.t initial value of IKK
            oed_a(k)=trace(inv(FIM));           %A-optimal criterion
            oed_d(k)=det(FIM);                  %D-optimal criterion
            oed_e(k)=min(lambda);               %E-optimal criterion
            oed_me(k)=max(lambda)/min(lambda);  %Modified E-optimal criterion         
  end
  
  
%################# PLOTTING ###############################
figure
    plot(IKK,oed_a,'r-');
    hold on;
    [a,b]=min(oed_a);
    IKK1=b/10*IKK_normal*ones(size(oed_a));
    plot(IKK1,oed_a,'b--');
    text(0.1,2e-3,['IKK=' num2str(b/10*IKK_normal) '\muM']);
    xlabel('Initial concentration of IKK (\muM)');
    ylabel('Trace(FIM^-^1)');
    legend('A-optimal criterion');
   
figure
    plot(IKK,oed_d,'r-');
    hold on;
    [a,b]=max(oed_d);
    IKK1=b/10*IKK_normal*ones(size(oed_d));
    plot(IKK1,oed_d,'b--');
    text(0.1,10e25,['IKK=' num2str(b/10*IKK_normal) '\muM']);
    xlabel('Initial concentration of IKK (\muM)');
    ylabel('Det(FIM)');
    legend('D-optimal criterion');
   
figure
    plot(IKK,oed_e,'r-');
    hold on;
    [a,b]=max(oed_e);
    IKK1=b/10*IKK_normal*ones(size(oed_e));
    plot(IKK1,oed_e,'b--');
    text(0.1,5e6,['IKK=' num2str(b/10*IKK_normal) '\muM']);
    xlabel('Initial concentration of IKK (\muM)');
    ylabel('\lambda_m_i_n(FIM)');
    legend('E-optimal criterion');
   
figure
    plot(IKK,oed_me,'r-');
    hold on;
    [a,b]=min(oed_me);
    IKK1=b/10*IKK_normal*ones(size(oed_me));
    plot(IKK1,oed_me,'b--');
    text(0.1,8e5,['IKK=' num2str(b/10*IKK_normal) '\muM']);
    xlabel('Initial concentration of IKK (\muM)');
    ylabel('\lambda_m_a_x(FIM)/\lambda_m_i_n(FIM)');
    legend('modified E-optimal criterion');
%######################################################################################################
%######################################################################################################





另外
parameter_hoff
%######################################################################################################
%######################################################################################################
num=600; %the number of sampling

%############## Simulation time points ##################################
t0=0; %time when IKK is being introduced into the system
tw=6000; %length of IKK stimulation
tau=(tw-t0)/num; %defines the period (in seconds) in between to successive data points

%############### Initial condition for total NF-kB ######################
y0=zeros(1,10); %Initial condition set to zero for each reaction participant
NF=0.9;
EF=1.5e-5;
SF=0.8;
y0(8)=NF; %NF-kB is given in cytoplasm
y0(1)=1.5e-5;
y0(6)=0.8;
Source=1;

%############### Parametrization for system of ODE's ####################
par=zeros(11,1);
par(1,1)=100000;
par(2,1)=1000;
par(3,1)=100;
par(4,1)=10000;
par(5,1)=50000;
par(6,1)=200;
par(7,1)=1000;
par(8,1)=20000;
par(9,1)=5000;
par(10,1)=100;
par(11,1)=500;

%######################################################################################################
%######################################################################################################


model_hoff
%######################################################################################################
function dy=model_hoff(t,y,par,Source)
dy=zeros(10,1);

par(1,1)   = 100000;
par(2,1)  = 1000;
par(3,1)   = 100;
par(4,1)  = 10000;
par(5,1)   = 50000;
par(6,1)  = 200;
par(7,1)   = 1000;
par(8,1)  = 20000;
par(9,1)  = 5000;
par(10,1)  = 100;
par(11,1)   = 500;

dy(1)=-par(1,1)*y(1)*y(6)+par(2,1)*y(2)+par(7,1)*y(4)-par(8,1)*y(1)*y(9)+par(11,1)*y(5);
dy(2)=par(1,1)*y(1)*y(6) - par(2,1)*y(2) - par(3,1)*y(2) + par(4,1)*y(3)*y(7);
dy(3)=par(3,1)*y(2) - par(4,1)*y(3)*y(7) - par(5,1)*y(3)*y(8) + par(6,1)*y(4) - par(9,1)*y(3) + par(10,1)*y(5);
dy(4)=par(5,1)*y(3)*y(8) - par(6,1)*y(4) - par(7,1)*y(4) + par(8,1)*y(1)*y(9);
dy(5)=par(9,1)*y(3) - par(10,1)*y(5) - par(11,1)*y(5);
dy(6)=-par(1,1)*y(1)*y(6) + par(2,1)*y(2);
dy(7)=par(3,1)*y(2) - par(4,1)*y(3)*y(7);
dy(8)=-par(5,1)*y(3)*y(8) + par(6,1)*y(4);
dy(9)=par(7,1)*y(4) - par(8,1)*y(1)*y(9);
dy(10)=par(11,1)*y(5);

%######################################################################################################
%######################################################################################################

怎么错了?  希望帮下忙
谢谢!!!
回复
分享到:

使用道具 举报

 楼主| 发表于 2011-3-20 19:51 | 显示全部楼层
回复 1 # lphtoby 的帖子

出现的错误是这样的
??? Attempted to access y(6); index out of bounds because numel(y)=1.

Error in ==> model_hoff at 17
dy(1)=-par(1,1)*y(1)*y(6)+par(2,1)*y(2)+par(7,1)*y(4)-par(8,1)*y(1)*y(9)+par(11,1)*y(5);

Error in ==> funfun\private\odearguments at 110
f0 = feval(ode,t0,y0,args{:});   % ODE15I sets args{1} to yp0.

Error in ==> ode15s at 227
[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, ...

Error in ==> IKK_optimal at 22
        [T,Y]=ode15s(@model_hoff,tspan,Y01,[],par,Source);

衷心感谢帮忙!!!
发表于 2011-3-20 20:20 | 显示全部楼层
本帖最后由 john152 于 2011-3-20 20:21 编辑

这问题提的
直接上来就是一大段程序
就算是让别人帮你改程序,你也得告知这是什么样的程序,是干嘛的吧?
发表于 2011-3-20 22:42 | 显示全部楼层
你可以试试这样:
function dy=model_hoff(t,y,par,Source)
dy=zeros(10,1)
y = zeros(10 , 1) ; % 加上去的!
我没试,你试试!
 楼主| 发表于 2011-3-20 23:28 | 显示全部楼层
回复 3 # john152 的帖子

这是优化设计的程序
 楼主| 发表于 2011-3-20 23:29 | 显示全部楼层
回复 3 # john152 的帖子

一共10个状态 11个参数 初始值已经知道  谢谢!!
 楼主| 发表于 2011-3-20 23:41 | 显示全部楼层
回复 4 # zhouyang664 的帖子

首先谢谢你的帮助!!
我加上
y=zero(10,1)
y(1)=1.5e-5;
y(2)=0;
y(3)=0;
y(4)=0;
y(5)=0;
y(6)=0.8;
y(7)=0;
y(8)=0.9;
y(9)=0;
y(10)=0;                %赋予初始值
结果运行提示为
??? Error using ==> funfun\private\odearguments
Solving MODEL_HOFF requires an initial condition vector of length 10.

Error in ==> ode15s at 227
[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, ...

Error in ==> IKK_optimal at 22
        [T,Y]=ode15s(@model_hoff,tspan,Y01,[],par,Source);

希望你有时间看看 行吗?
在此感谢您的帮忙!!!
发表于 2011-3-21 20:21 | 显示全部楼层
试了试,没整出来,没有太多时间,还是楼主自己好好整整吧!
发表于 2011-8-2 19:26 | 显示全部楼层
楼主现在整明白了吗?问题解决了没?
发表于 2011-8-3 00:50 | 显示全部楼层
??? Attempted to access y(6); index out of bounds because numel(y)=1.

报错都已经告知y仅为纯量, 那来的y(6), y(5) ...
发表于 2011-9-13 11:22 | 显示全部楼层
你好  我跟你的问题机会一样  请问你的解决了吗

点评

http://www.chinavib.com/thread-105682-1-1.html  发表于 2011-9-14 10:13
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-14 14:49 , Processed in 0.062127 second(s), 18 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

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