Plotting the mean and standard derivation of many stair-step graphs
3 views (last 30 days)
Show older comments
I want to find the mean and standard derivation of several stairstep graphs with different length of steps.
Here is my function with some random variables.
clear all
%set parameter values
am=1;
bm=1;
ap=1;
bp=1;
ohm=100;
%set number of reactions
Tend=4000;
%set initial conditions for # molecules and for time t
X=zeros(Tend,2);
t=zeros(Tend,1);
X(1,:)=[0 0]; t(1)=0;
j=1;
while t(j)<10
%calculate propensities
p1=am*ohm;
p2=bm*X(j,1);
p3=ap*X(j,1);
p4=bp*X(j,2);
psum=p1+p2+p3+p4;
%select next reaction
mu=rand(1);
z1=0;z2=0;z3=0;z4=0;
if 0 <= mu && mu < p1/psum % M generated
z1=1;
else if p1/psum <= mu && mu < (p1+p2)/psum % M degraded
z2=1;
else if (p1+p2)/psum <= mu && mu < (p1+p2+p3)/psum % P generated
z3=1;
else % P degraded
z4=1;
end
end
end
%update state
X(j+1,1) = X(j,1) + z1 -z2;
X(j+1,2) = X(j,2) + z3 -z4;
%update time
t(j+1) = t(j) + log(1/rand(1))/psum;
%update counter
j=j+1;
end
subplot(2,1,1)
stairs(t(1:j), X(1:j,1), 'k', 'linewidth', 2)
xlim([0 10])
ylabel('number of molecules')
xlabel('time')
title('M')
subplot(2,1,2)
stairs(t(1:j), X(1:j,2), 'g', 'linewidth', 2)
xlim([0 10])
ylabel('number of molecules')
xlabel('time')
title('P')
This is the result of the first trial, and
this is the one of the second trial.
You can see that as the graph differs, j and the length of steps also differs everytime the function is plotted.
So I cannot use stdshade function since it assumes that 'the observation' occurs with uniform steps of the time.
I tried to find the value of M and P at t=0.01*n to make a new matrix, and then use stdshade.
while k<=1000
while t(j)<10
%calculate propensities
p1=am*ohm;
p2=bm*X(j,1);
p3=ap*X(j,1);
p4=bp*X(j,2);
psum=p1+p2+p3+p4;
%select next reaction
mu=rand(1);
z1=0;z2=0;z3=0;z4=0;
if 0 <= mu && mu < p1/psum % M generated
z1=1;
else if p1/psum <= mu && mu < (p1+p2)/psum % M degraded
z2=1;
else if (p1+p2)/psum <= mu && mu < (p1+p2+p3)/psum % P generated
z3=1;
else % P degraded
z4=1;
end
end
end
%update state
X(j+1,1) = X(j,1) + z1 -z2;
X(j+1,2) = X(j,2) + z3 -z4;
%update time
t(j+1) = t(j) + log(1/rand(1))/psum;
%update counter
j=j+1;
end
f1=stairs(t(1:j),X(1:j,1));
f2=stairs(t(1:j),X(1:j,2));
XData1=get(get(gca,'f1'),'XData');
YData1=get(get(gca,'f1'),'YData');
XData2=get(get(gca,'f2'),'XData');
YData2=get(get(gca,'f2'),'YData');
while n<=1000
M1(k,n)=YData1(XData1==0.01*n);
M2(k,n)=YData2(XData2==0.01*n);
n=n+1;
end
n=1;
j=1;
end
stdshade(M1,0:0.01:10)
stdshade(M2,0:0.01:10)
But it shows the error at 'matlab.graphics.axis.Axes/get'.
The expected result is as follows:
Can you help me? Thanks.
2 Comments
dpb
on 14 Nov 2021
We can't run your function nor did you even show us what it produces -- not a lot to be done with such limited info.
That said, why can't you just use interp1 to create a set of uniform points for all the curves and then average those?
Answers (0)
See Also
Categories
Find more on 2-D and 3-D Plots in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!