How can we create kappa and delta in the next model using matlab?

Question

I have the following stochastic model that describes the evolution of a process (Y) in space and time. Ds and Dt are the area in space (2D with the x and y axes) and time (1D with the t axis). This model is commonly known as a mixed effect model or variation model.

I am currently developing Y as follows:

%# Time parameters T=1:1:20; % input nT=numel(T); %# Grid and model parameters nRow=100; nCol=100; [Grid.Nx,Grid.Ny,Grid.Nt] = meshgrid(1:1:nCol,1:1:nRow,T); xPower=0.1; tPower=1; noisePower=1; detConstant=1; deterministic_mu = detConstant.*(((Grid.Nt).^tPower)./((Grid.Nx).^xPower)); beta_s = randn(nRow,nCol); % mean-zero random effect representing location specific variability common to all times gammaTemp = randn(nT,1); for t = 1:nT gamma_t(:,:,t) = repmat(gammaTemp(t),nRow,nCol); % mean-zero random effect representing time specific variability common to all locations end var=0.1;% noise has variance = 0.1 for t=1:nT kappa_st(:,:,t) = sqrt(var)*randn(nRow,nCol); end for t=1:nT Y(:,:,t) = deterministic_mu(:,:,t) + beta_s + gamma_t(:,:,t) + kappa_st(:,:,t); end 

My questions:

• How to produce a delta in the expression for Y and the difference in kappa and delta?
• Help explain using some illustrations using Matlab if I create Y correctly?

Please let me know if you need more information / explanation. Thanks.

Answer 1

Your samples of beta, gamma and kappa implementations, as if they are white (for example, their values ​​for each (x, y, t) are independent). Descriptions of terms suggest that this is not so. It seems that the delta should capture white noise, while other terms capture correlations in their respective regions. for example, there is a non-zero correlation between gamma (t_1) and gamma (t_1 + 1).

If you want to model gamma as a stationary Gaussian Markov process with var_g dispersion and cor_g correlation between gamma (t) and gamma (t + 1), you can use something like

 gamma_t = nan( nT, 1 ); gamma_t(1) = sqrt(var_g)*randn(); K_g = cor_g/var_g; K_w = sqrt( (1-K_g^2)*var_g ); for t = 2:nT, gamma_t(t) = K_g*gamma_t(t-1) + K_w*randn(); end gamma_t = reshape( gamma_t, [ 1 1 nT ] ); 

The formulas that I used for the K_g and K_w gains in the above code (and the initialization gamma_t (1)) create the desired stationary variance \ sigma ^ 2_0 and one-step covariance \ sigma ^ 2_1:

Note that the implementation above assumes that you will later summarize the terms with bsxfun to do a “repmat” for you:

 Y = bsxfun( @plus, deterministic_mu + kappa_st + delta_st, beta_s ); Y = bsxfun( @plus, Y, gamma_t ); 

Please note that I have not tested the above code, so you must confirm with a sample that it really produces zero noise with a given dispersion and covariance between adjacent samples. To try the beta version, the same procedure can be extended to two dimensions, but the principles are essentially the same. I suspect that kappa should be similarly modeled as a Markov Gaussian process, but in all three dimensions and with lower dispersion to represent higher order effects not captured in mu, beta and gamma.

The delta should be zero mean stationary white noise. Assuming this is Gaussian with noisePower dispersion, you could try it out with

 delta_st = sqrt(noisePower)*randn( [ nRows nCols nT ] ); 

Answer 2

First, I rewrote your code to make it more efficient. I see that you are creating linear exploded meshes for x , y and t and doing calculations for all points in that grid. This approach has severe limitations on the maximum achievable grid resolution, since the 3D grid (and all variables defined with it) can consume a very large amount of memory if the resolution increases. If the model you are implementing will increase in complexity and size (this often happens), I would suggest that you translate all this into a function that takes matrix / vector inputs for s and t , which will be a little more flexible in this regard - processing " blocks of data that would otherwise not fit into memory would be much simpler.

Then I generated a delta_st member with rand instead of randn , since the noise should be white. Now I’m very unsure of this last one, and I didn’t have time to read the article that you contacted - can you tell me on which pages I can find the relevant sections for delta_st ?

Now the code:

 %# Time parameters T = 1:1:20; % input nT = numel(T); %# Grid and model parameters nRow = 100; nCol = 100; % noise has variance = 0.1 var = 0.1; xPower = 0.1; tPower = 1; noisePower = 1; detConstant = 1; [Grid.Nx,Grid.Ny,Grid.Nt] = meshgrid(1:nCol,1:nRow,T); % deterministic mean deterministic_mu = detConstant .* Grid.Nt.^tPower ./ Grid.Nx.^xPower; % mean-zero random effect representing location specific % variability common to all times beta_s = repmat(randn(nRow,nCol), [1 1 nT]); % mean-zero random effect representing time specific % variability common to all locations gamma_t = bsxfun(@times, ones(nRow,nCol,nT), randn(1, 1, nT)); % mean zero random effect capturing the spatio-temporal % interaction not found in the larger-scale deterministic mu kappa_st = sqrt(var)*randn(nRow,nCol,nT); % mean zero random effect representing the micro-scale % spatio-temporal variability that is modelled by white % noise (iid at different time steps) in Ds·Dt delta_st = noisePower * (rand(nRow,nCol,nT)-0.5); % Final result: Y = deterministic_mu + beta_s + gamma_t + kappa_st + delta_st;