论文求助：关于用Matlab实现分段函数求变上限积分的问题

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

大家好，我正在写论文，遇到了如下的问题，很棘手，比较焦急，请大家指教！

• syms f(t) t m；
• thetatL(i-1)=-0.5;thetatH(i-1)=1;
• xt(i)=0.5;
• mt(j,i)=1;
• f(t)=feval(symengine,'piecewise','[t<thetatL(i-1), 0]','[t>=thetatL(i-1)&t<=thetatH(i-1), (normcdf(xt(i)*(t-mt(j,i)))-normcdf(xt(i)*(thetatL(i-1)-mt(j,i))))/(normcdf(xt(i)*(thetatH(i-1)-mt(j,i)))-normcdf(xt(i)*(thetatL(i-1)-mt(j,i))))]','[t>thetatH(i-1),1]') ；
• m=int(f(t),-Inf,x);
• m=matlabFunction(m);
• eq=@(x) m-1/2;
• sol=fsolve(eq,0);
• xcrijt(i)=sol;
  

[color=rgb(51, 102, 153) !important]复制代码


这段之前是有一个关于i和j的循环，在上一次中算出了xt(i)、mt(j,i)、thetatL(i-1) 和 thetatH(i-1) 各自的值，为了简便我在这里假设了它们的值。

f(t)其实是一个条件分布函数，依据 i-1 次的thetatL(i-1) 和 thetatH(i-1)来确定 t的分布，为了便于大家看清，我把它的形式附在下图里：



我的目标是在分布更新后求出满足Prob{t<=x}=1/2的x的值，也就是代码里的xcrijt(i)。


运行完代码后，我得到下面一堆错误。。。
警告: Cannot check whether the integrand is defined everywhere on the integration interval. [defint_b_gr_a]

错误使用 mupadmex
Error in MuPAD command: Cannot generate code for object 'thetatL(- 1 + I) in R_'. [_in::CF]

出错 sym/matlabFunction>mup2mat (line 316)
res = mupadmex('symobj::generateMATLAB',r.s,0);

出错 sym/matlabFunction>mup2matcell (line 304)
    r = mup2mat(c{1});

出错 sym/matlabFunction (line 123)
    body = mup2matcell(funs);

出错 Untitled4 (line 85)
   m=matlabFunction(m);


新手入门，请大家帮帮忙，谢谢大家的指教！感激不尽！

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享0 收藏1 回帖

关键词：matlab实现 MATLAB matla atlab 论文求助 分段函数 变限积分 matlab

function pv = ft_ex
th_L=-0.5; th_H=1; xt = 0.5; mt = 1;
PL = normcdf(xt*(th_L-mt)); PH = normcdf(xt*(th_H-mt));
ft = @(t) 0+(t>=th_L & t<=th_H).*(normcdf(xt*(t-mt))-PL)/(PH-PL)+(t>th_H);
Target = @(x) integral(ft, -inf, x)- 0.5;
pv = fzero(@(p) Target(p), [th_L,th_H]);
.......
pv=0.7642

jyliao 发表于 2016-8-21 08:56
function pv = ft_ex
th_L=-0.5; th_H=1; xt = 0.5; mt = 1;
PL = normcdf(xt*(th_L-mt)); PH = normcdf( ...

您好！谢谢您的回复！我发现之前的做法有错误，这两天受您的启发改了改，但是仍然有错误，能帮我看一下吗？
%前面由循环语句得出来ft的形式是

ft = @(t) 0+(t>th_H).*(normcdf((t-mt(j,i))/xt(i))-PH)/(1-PH); %或者
ft = @(t) (t>=th_L)+(t<th_L).*(normcdf((t-mt(j,i)))/xt(i))/PL; %或者
ft = @(t) 0+(t>=th_L & t<=th_H).*(normcdf((t-mt(j,i))/xt(i))-PL)/(PH-PL)+(t>th_H);


gt=@(t,h) ft(t)*(normpdf(sqrt(ys(i))*(h-t)));

Tt = @(h) sqrt(ys(i))*integral(@(t) gt, -inf, inf)- 0.5; %我想就gt这个二元函数关于t取期望，得到一个h的函数，带到下面的方程里面找零点

pv = fzero(@(p) Tt(p), -10000);

%但是报错了：

错误使用 fzero (line 289)
FZERO 无法继续，因为用户提供的 function_handle ==> @(p)0.5-Tt(p) 失败，出现下面的错误。


输入函数必须返回 'double' 或 'single' 值。找到 'function_handle'。

出错 trial2 (line 131)
    pv = fzero(@(p) 0.5-Tt(p), -10000);

Tt = @(h) sqrt(ys(i))*integral(@(t) gt, -inf, inf)- 0.5;
h is not in the integral expression

jyliao 发表于 2016-8-24 20:24
Tt = @(h) sqrt(ys(i))*integral(@(t) gt, -inf, inf)- 0.5;
h is not in the integral expression

后来我放弃了这种求法，转成t=-1000:0.005:1000，作两个pdf的conv了，问题解决了但是感觉计算量变得好大，加上程序的其它循环笔记本有点吃不消了，i5 6300u跑一个半小时才跑出来三个数。。。