%% Process equation x[k] = sys(k, x[k-1], u[k]);% untuk state vector
nx = 10; % number of states
sys=@(k, xkm1, uk) xkm1./2 + 25.*xkm1/(1+xkm1.^2) + 8*cos(1.2*k) + uk; % (returns column vector)
%% Observation equation y[k] = obs(k, x[k], v[k]);
ny = 1; % number of observations
obs = @(k, xk,vk) 100 + vk; % (returns column vector)
%% PDF of process noise and noise generator function
nu = 10; % size of the vector of process noise
sigma_u = sqrt(0.0000001);
p_sys_noise = @(u) normpdf(u, 0, sigma_u);
gen_sys_noise = @(u) normrnd(0, sigma_u); % sample from p_sys_noise (returns column vector)
Notice that gen_sys_noise accepts a parameter but completely ignores it -- not, for example, using it to determine the size of the desired output.
%% PDF of observation noise and noise generator function
nv = 1; % size of the vector of observation noise
sigma_v = sqrt(0.0000001);
p_obs_noise = @(v) normpdf(v, 0, sigma_v);
gen_obs_noise = @(v) normrnd(0, sigma_v); % sample from p_obs_noise (returns column vector)
Likewise gen_obs_noise ignores its parameter.
%% Initial PDF
% p_x0 = @(x) normpdf(x, 0,sqrt(10)); % initial pdf
gen_x0 = @(x) [95;96;97;98;94;93;97;99;98;100]+ones(10,1)*normrnd(0, sqrt(0.0000001)); % sample from p_x0 (returns column vector)
%% Transition prior PDF p(x[k] | x[k-1])
% (under the suposition of additive process noise)
% p_xk_given_xkm1 = @(k, xk, xkm1) p_sys_noise(xk - sys(k, xkm1, 0));
%% Observation likelihood PDF p(y[k] | x[k])
% (under the suposition of additive process noise)
p_yk_given_xk = @(k, yk, xk) p_obs_noise(yk - obs(k, xk, 0));
%% %% Number of time steps
T = 100;
%% Separate memory space
x = zeros(nx,T); y = zeros(ny,T);
u = zeros(nu,T); v = zeros(nv,T);
%% Simulate system
xh0 = 0; % initial state
u(:,1) = 0; % initial process noise
v(:,1) = gen_obs_noise(sigma_v); % initial observation noise
x(:,1) = xh0;
y(:,1) = obs(1, xh0, v(:,1));
for k = 2:T
% here we are basically sampling from p_xk_given_xkm1 and from p_yk_given_xk
u_temp = gen_sys_noise(); % simulate process noise
v_temp = gen_obs_noise(); % simulate observation noise
x_temp = sys(k, x(:,k-1), u(:,k)); % simulate state
y_temp = obs(k, x(:,k), v(:,k)); % simulate observation
whos u v x y u_temp v_temp x_temp y_temp
disp('about to u_temp')
u(:,k) = u_temp;
disp('about to v_temp')
v(:,k) = v_temp;
disp('about to x_temp')
x(:,k) = x_temp;
disp('about to y_temp')
y(:,k) = y_temp;
disp('got through temp')
end
the gen_*_noise functions each return a scalar. You assign the scalar to a (:,k) of an array, and that is fine as far as MATLAB cares: it copies the scalar to each of the output locations. Whether that makes sense mathematically for your case is a different matter -- should you perhaps be using a different noise value for each row?
You are calling sys() with column vectors for your second and third parameter. I explained in detail to you in https://www.mathworks.com/matlabcentral/answers/1848368-meaning-of-sys-k-xkm1-uk-xkm1-2-25-xkm1-1-xkm1-2-8-cos-1-2-k-uk#answer_1096918 that the / operator is the matrix-right-divide operation
"A matrix matrix-right-divide another matrix of the same size is well defined, and will be roughly the same as (25*xkm1) * pinv(1+xkm1^2) where pinv is pseudo-inverse and * is the inner-product operation. You will get out a matrix the same size as xkm1 ... provided that xkm1 was a square matrix"
but your xkm1 is not a square matrix. When you divide a row vector by another row vector the same size, the result is a square matrix that has as many rows and columns as there are rows in the row vectors.
What you should have understood from my long explaination is that you should not be using the MATLAB / operator, except possibly in cases where the right hand side is a scalar. A/2 is not going to cause you problems, but A/B where B is non-scalar is seldom what you want to do.
Only use A/B for cases where B is guaranteed scalar, or where you mean something algebraically equivalent to innerproduct between A and matrix inverse of B.
