From 0a9438ee079d2a8d4446e664f64abf680bb08afa Mon Sep 17 00:00:00 2001 From: Arno Solin Date: Sat, 4 May 2019 22:26:46 +0300 Subject: [PATCH] Book example codes. --- LICENSE | 4 +- README.md | 179 +++++++- matlab/ch02_ex09_numerical_solution_of_odes.m | 108 +++++ matlab/ch03_ex07_time_varying_oscillator.m | 84 ++++ matlab/ch03_ex10_stochastic_spring_model.m | 101 +++++ .../ch03_fig02_two_views_of_brownian_motion.m | 76 ++++ matlab/ch03_fig10_white_noise.m | 24 ++ matlab/ch04_ex05_ou_process.m | 58 +++ matlab/ch04_fig01_brownian_motion.m | 22 + matlab/ch07_ex03_karhunen_loeve_series.m | 48 +++ matlab/ch07_ex12_doobs_h_transform.m | 60 +++ matlab/ch07_ex17_feynman_kac.m | 309 ++++++++++++++ matlab/ch08_ex06_weak_itotaylor.m | 169 ++++++++ matlab/ch08_ex11_ode_solvers.m | 113 +++++ matlab/ch08_ex15_duffing_van_der_pol.m | 170 ++++++++ matlab/ch08_ex21_leapfrog_verlet.m | 181 ++++++++ matlab/ch08_ex24_exact_simulation.m | 134 ++++++ matlab/ch09_ex06_linearizations_for_benes.m | 398 ++++++++++++++++++ matlab/ch09_ex07_gaussian_approximation.m | 83 ++++ matlab/ch09_ex18_hermite_expansion.m | 67 +++ matlab/ch09_ex20_discretized_fpk_for_benes.m | 109 +++++ matlab/ch09_ex23_pathwise_series_expansion.m | 132 ++++++ matlab/ch10_ex17_benes_daum_ekf_erts.m | 299 +++++++++++++ matlab/ch10_ex19_ou_filtering_smoothing.m | 269 ++++++++++++ matlab/ch11_ex05_ou_parameter_estimation.m | 381 +++++++++++++++++ matlab/ch12_ex06_gp_regression.m | 226 ++++++++++ matlab/ch12_ex12_gp_approximation_of_drift.m | 111 +++++ matlab/euler.m | 49 +++ matlab/eulermaruyama.m | 63 +++ matlab/eulermaruyama_weak.m | 74 ++++ matlab/herm_exp.m | 53 +++ matlab/heun.m | 52 +++ matlab/impliciteuler.m | 51 +++ matlab/leapfrog.m | 73 ++++ matlab/lti_disc.m | 67 +++ matlab/ndhist.m | 173 ++++++++ matlab/newquiver.m | 123 ++++++ matlab/rejectionsample.m | 54 +++ matlab/rk4simple.m | 56 +++ matlab/srkS10scalarnoise.m | 81 ++++ matlab/srkW20.m | 111 +++++ matlab/w20scalar.m | 97 +++++ 42 files changed, 5088 insertions(+), 4 deletions(-) create mode 100644 matlab/ch02_ex09_numerical_solution_of_odes.m create mode 100644 matlab/ch03_ex07_time_varying_oscillator.m create mode 100644 matlab/ch03_ex10_stochastic_spring_model.m create mode 100644 matlab/ch03_fig02_two_views_of_brownian_motion.m create mode 100644 matlab/ch03_fig10_white_noise.m create mode 100644 matlab/ch04_ex05_ou_process.m create mode 100644 matlab/ch04_fig01_brownian_motion.m create mode 100644 matlab/ch07_ex03_karhunen_loeve_series.m create mode 100644 matlab/ch07_ex12_doobs_h_transform.m create mode 100644 matlab/ch07_ex17_feynman_kac.m create mode 100644 matlab/ch08_ex06_weak_itotaylor.m create mode 100644 matlab/ch08_ex11_ode_solvers.m create mode 100644 matlab/ch08_ex15_duffing_van_der_pol.m create mode 100644 matlab/ch08_ex21_leapfrog_verlet.m create mode 100644 matlab/ch08_ex24_exact_simulation.m create mode 100644 matlab/ch09_ex06_linearizations_for_benes.m create mode 100644 matlab/ch09_ex07_gaussian_approximation.m create mode 100644 matlab/ch09_ex18_hermite_expansion.m create mode 100644 matlab/ch09_ex20_discretized_fpk_for_benes.m create mode 100644 matlab/ch09_ex23_pathwise_series_expansion.m create mode 100644 matlab/ch10_ex17_benes_daum_ekf_erts.m create mode 100644 matlab/ch10_ex19_ou_filtering_smoothing.m create mode 100644 matlab/ch11_ex05_ou_parameter_estimation.m create mode 100644 matlab/ch12_ex06_gp_regression.m create mode 100644 matlab/ch12_ex12_gp_approximation_of_drift.m create mode 100644 matlab/euler.m create mode 100644 matlab/eulermaruyama.m create mode 100644 matlab/eulermaruyama_weak.m create mode 100644 matlab/herm_exp.m create mode 100644 matlab/heun.m create mode 100644 matlab/impliciteuler.m create mode 100644 matlab/leapfrog.m create mode 100644 matlab/lti_disc.m create mode 100644 matlab/ndhist.m create mode 100644 matlab/newquiver.m create mode 100644 matlab/rejectionsample.m create mode 100644 matlab/rk4simple.m create mode 100644 matlab/srkS10scalarnoise.m create mode 100644 matlab/srkW20.m create mode 100644 matlab/w20scalar.m diff --git a/LICENSE b/LICENSE index 408246c..25c1505 100644 --- a/LICENSE +++ b/LICENSE @@ -1,6 +1,6 @@ MIT License -Copyright (c) 2019 AaltoML +Copyright (c) 2019 Simo Särkkä and Arno Solin Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal @@ -18,4 +18,4 @@ FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -SOFTWARE. +SOFTWARE. \ No newline at end of file diff --git a/README.md b/README.md index fd522ba..71bf6fc 100644 --- a/README.md +++ b/README.md @@ -1,2 +1,177 @@ -# SDE -Example codes for the book Applied Stochastic Differential Equations +# Applied Stochastic Differential Equations + +[Simo Särkkä](https://users.aalto.fi/~ssarkka/) · [Arno Solin](http://arno.solin.fi) + +Example codes for the book: + +* Simo Särkkä and Arno Solin (2019). **Applied Stochastic Differential Equations**. Cambridge University Press. Cambridge, UK. + +## Summary + +The book *Applied Stochastic Differential Equations* gives a gentle introduction to stochastic differential equations (SDEs). The low learning curve only assumes prior knowledge of ordinary differential equations and basic concepts of statistic, together with understanding of linear algebra, vector calculus, and Bayesian inference. The book is mainly intended for advanced undergraduate and graduate students in applied mathematics, signal processing, control engineering, statistics, and computer science. + +The worked examples and numerical simulation studies in each chapter illustrate how the theory works in practice and can be implemented for solving the problems. To promote hands-on work with the methods, we provide the MATLAB􏰀 source code for reproducing the example results in the book. The code examples have been grouped by chapter, and some pointers to example and figure numbers in the book are given below. + +## Codes for the examples + +All codes for the examples (excluding pen and paper examples which do not have any code associated to them) and figures (those requiring numerical simulation) in the book are provided below. Some entries cover multiple examples/figures in the same chapter, in the case of which the code file naming follows the first example. + +### Chapter 2: Numerical solution of ODEs + +This experiment replicates the results in Example 2.9 (Numerical solution of ODEs) in the book (Fig. 2.1). + +```ch02_ex09_numerical_solution_of_odes.m``` + +### Chapter 3: Two views of Brownian motion + +This experiment replicates the results in Figure 3.2 in the book. + +```ch03_fig02_two_views_of_brownian_motion.m``` + +### Chapter 3: Time-varying oscillator + +This experiment replicates the results in Example 3.7 (Heart and breathing tracking in the brain) in the book (Fig. 3.8). + +```ch03_ex07_time_varying_oscillator.m``` + +### Chapter 3: Stochastic spring model + +This experiment replicates the results in Example 3.10 (Stochastic spring model) in the book (Figs. 3.9 and 3.11). + +```ch03_ex10_stochastic_spring_model.m``` + +### Chapter 3: White noise + +This experiment replicates the results in Figure 3.10 in the book. + +```ch03_fig10_white_noise.m``` + +### Chapter 4: Brownian motion + +This experiment replicates the results in Figure 4.1 in the book. + +```ch04_fig01_brownian_motion.m``` + +### Chapter 4: Solution of the Ornstein–Uhlenbeck process + +This experiment replicates the results in Example 4.5 (Solution of the Ornstein–Uhlenbeck process) in the book (Fig. 4.2). + +```ch04_ex05_ou_process.m``` + +### Chapter 7: Karhunen–Loeve series of Brownian motion + +This experiment replicates the results in Example 7.3 (Karhunen–Loeve series of Brownian motion) in the book (Fig. 7.1). + +```ch07_ex03_karhunen_loeve_series.m``` + +### Chapter 7: Doob's h-transform + +This experiment replicates the results in Example 7.12 (Conditioned Ornstein–Uhlenbeck process) in the book (Fig. 7.2). + +```ch07_ex12_doobs_h_transform.m``` + +### Chapter 7: Feynman-Kac formula + +This experiment replicates the results in Example 7.17 (Solution of an elliptic PDE using SDE simulation) in the book (Fig. 7.3). + +```ch07_ex17_feynman_kac.m``` + +### Chapter 8: Weak Gaussian vs. weak three-point approximations + +This experiment replicates the results in Example 8.6 (Simulating from a trigonometric nonlinear SDE) in the book (Fig. 8.1). + +```ch08_ex06_weak_itotaylor.m``` + +### Chapter 8: Comparison of Runge–Kutta schemes + +This experiment replicates the results in Example 8.11 (Comparison of ODE solvers) in the book (Fig. 8.2). + +```ch08_ex11_ode_solvers.m``` + +### Chapter 8: Duffing van der Pol model + +This experiment replicates the results in Example 8.15 (Duffing van der Pol oscillator) in the book (Figs. 8.3–8.5). + +```ch08_ex15_duffing_van_der_pol.m``` + +### Chapter 8: Leapfrog Verlet + +This experiment replicates the results in Example 8.21 (Leapfrog solution to the spring model) in the book (Fig. 8.6). + +```ch08_ex21_leapfrog_verlet.m``` + +### Chapter 8: Exact simulation of sine diffusion + +This experiment replicates the results in Example 8.24 (Exact simulation of sine diffusion) in the book (Fig. 8.7). + +```ch08_ex24_exact_simulation.m``` + +### Chapter 9: Linearizations and approximations for the Beneš model + +This experiment replicates the results in Example 9.6 (Linearization and Gauss–Hermite approximations, shown in Fig. 9.1) and Example 9.14 (Local linearization vs. Gaussian approximations, shown in Fig. 9.3). + +```ch09_ex06_linearizations_for_benes.m``` + +### Chapter 9: Gaussian approximation + +This experiment replicates the results in Example 9.7 (Gaussian approximation of a nonlinear trigonometric SDE) in the book (Fig. 9.2). + +```ch09_ex07_gaussian_approximation.m``` + +### Chapter 9: Hermite expansion of Beneš SDE + +This experiment replicates the results in Example 9.18 (Hermite expansion of Beneš SDE) in the book (Fig. 9.4). + +```ch09_ex18_hermite_expansion.m``` + +### Chapter 9: Discretized FPK for the Beneš SDE + +This experiment replicates the results in Example 9.20 (Discretized FPK for the Beneš SDE) in the book (Fig. 9.5). + +```ch09_ex20_discretized_fpk_for_benes.m``` + +### Chapter 9: Pathwise series expansion of the Beneš SDE + +This experiment replicates the results in Example 9.23 (Pathwise series expansion of Beneš SDE) in the book (Fig. 9.6). + +```ch09_ex23_pathwise_series_expansion.m``` + +### Chapter 10: Beneš-Daum and EKF/ERTS examples + +This experiment replicates the results in Example 10.17 (Beneš–Daum filter, Fig. 10.2), Example 10.26 (Continuous-discrete EKF solution to the Beneš–Daum problem, Fig. 10.4), and Example 10.38 (Extended RTS solution to Beneš and Beneš–Daum filtering problems, Fig. 10.6) in the book. + +```ch10_ex17_benes_daum_ekf_erts.m``` + +### Chapter 10: Ornstein-Uhlenbeck filtering and smoothing + +This experiment replicates the results in Example 10.19 (Kalman filter for the Ornstein–Uhlenbeck model) and Example 10.21 (Continuous-discrete Kalman filter for the Ornstein–Uhlenbeck model), the results of which are shared in Figure 10.3. Additionally, the experiment also covers Example 10.29 (RTS smoother for the Ornstein–Uhlenbeck model) and Example 10.33 (Continuous-discrete/continuous RTS smoother for the Ornstein–Uhlenbeck model), the results of which are in Figure 10.5. + +```ch10_ex19_ou_filtering_smoothing.m``` + +### Chapter 11: Parameter estimation in Ornstein-Uhlenbeck + +This experiment replicates the results in Example 11.5 (Exact parameter estimation in the Ornstein–Uhlenbeck model, Fig. 11.1) and Example 11.9 (Approximate parameter estimation in the Ornstein–Uhlenbeck model, Fig. 11.2) in the book. + +```ch11_ex05_ou_parameter_estimation.m``` + +### Chapter 12: Batch and sequential solution to GP regression + +This experiment replicates the results in Example 12.6 (Batch GP regression, Fig. 12.1) and Example 12.11 (Sequential solution to GP regression, Fig. 12.2) in the book. + +```ch12_ex06_gp_regression.m``` + +### Chapter 12: GP approximation of the drift function (double-well) + +This experiment replicates the results in Example 12.12 (GP approximation of the double-well model) in the book (Fig. 12.3). + +```ch12_ex12_gp_approximation_of_drift.m``` + +## How to run? + +These codes have been tested under *Mathworks MATLAB R2018b* and *GNU Octave 4.4*. The code is pseudo-code like and tries to follow the presentation in the book. Thus they should be applicable for porting to other languages as well (consider this an exercise). + +## License + +Copyright Simo Särkkä and Arno Solin. + +This software is provided under the MIT License. See the accompanying LICENSE file for details. \ No newline at end of file diff --git a/matlab/ch02_ex09_numerical_solution_of_odes.m b/matlab/ch02_ex09_numerical_solution_of_odes.m new file mode 100644 index 0000000..1af892b --- /dev/null +++ b/matlab/ch02_ex09_numerical_solution_of_odes.m @@ -0,0 +1,108 @@ +%% Example 2.9: Numerical solution of ODEs +% +% Copyright: +% 2018 - Simo Särkkä and Arno Solin +% +% License: +% This software is provided under the MIT License. See the accompanying +% LICENSE file for details. + +%% Simulate model + + % Specify method (add your method here) + method_name = {'Euler','Heun','RK4'}; + method_handle = {@euler,@heun,@rk4simple}; + lineStyles = {'-','--','-.'}; + + % Set parameters + t1 = 10; + dt = .1; + g = 1; + v = 2; + q = 0.02; + x0 = [1;0]; + + % Dynamics + F = [0 1; -v^2 -g]; + f = @(x,t) F*x; + + % Exact solution + t_exact = linspace(0,t1,500); + x_exact = zeros(size(F,1),numel(t_exact)); + for j=1:numel(t_exact) + x_exact(:,j) = expm(F*t_exact(j))*x0; + end + + % Allocate space for approximate results + t = cell(length(method_name),1); + x = cell(length(method_name),1); + + % Run each method + for j=1:length(method_name) + [x{j},t{j}] = feval(method_handle{j},f,[0:dt:t1]',x0); + end + + % Show result + figure(1); clf; hold on + + % Plot exact + plot(t_exact,x_exact(1,:),'-','LineWidth',2,'Color',[.5 .5 .5]) + + % Plot numerical results + for j=1:length(method_name) + plot(t{j},x{j}(1,:),lineStyles{j},'Color','k') + end + + % Pimp up the plot + legend(['Exact' method_name]) + box on + xlabel('Time, $t$') + ylabel('$x_1$') + ylim([-.8 1]) + + + +%% Error plots + + % Step sizes + dt = logspace(-3,-1,8); + + % Errors + err = zeros(length(method_name),numel(dt)); + + % For each step length + for i=1:numel(dt) + + % Time steps + t = 0:dt(i):t1; + + % Exact result + x_exact = zeros(size(F,1),numel(t)); + for j=1:numel(t) + x_exact(:,j) = expm(F*t(j))*x0; + end + + % Run each method + for j=1:length(method_name) + + % Evaluate + x = feval(method_handle{j},f,t,x0); + + % Calcualte absolute error + err(j,i) = mean(abs(x(:)-x_exact(:))); + + end + end + + figure(2); clf; hold on + for j=1:length(method_name) + plot(dt,err(j,:),lineStyles{j},'Color','k') + end + plot(dt,err,'x','Color','k') + set(gca,'XScale','log','YScale','log') + box on + xlabel('$\Delta t$') + ylabel('Absolute error, $|\hat{\vx} - \vx|$') + set(gca,'YTick',10.^[-12 -8 -4 0], ... + 'YTickLabel',{'$10^{-12}$','$10^{-8}$','$10^{-4}$','$10^{0}$'}) + \ No newline at end of file diff --git a/matlab/ch03_ex07_time_varying_oscillator.m b/matlab/ch03_ex07_time_varying_oscillator.m new file mode 100644 index 0000000..87a4e20 --- /dev/null +++ b/matlab/ch03_ex07_time_varying_oscillator.m @@ -0,0 +1,84 @@ +%% Example 3.7: Stochastic resonators with time-varying frequency +% +% Copyright: +% 2018 - Simo Särkkä and Arno Solin +% +% License: +% This software is provided under the MIT License. See the accompanying +% LICENSE file for details. + +%% + + % Lock random seed + if exist('rng') % Octave doesn't have rng + rng(0,'twister') + else + randn('state',1); + end + + % Time + t = linspace(0,10,512); + + % Define the sinc function + sinc = @(x) sin(x)./x; + + % Set time-varying frequency trajectory + f = 1+.5*abs(t-5); + + % Parameters + H = [1 0]; + + figure(1); clf; hold on + + for j=1:4 + + % Allocate space + x = zeros(2,numel(t)); + z = x; + + % For harmonics + for i=1:2 + + % Diffusion + Qc = 10^-i; + + % Initial state + x(:,1) = [0; 1]; + + % Loop through time points + for k=2:numel(t) + + % Time delta + dt = t(k)-t(k-1); + + % Define F and L + F = [0 2*pi*i*f(k); -2*pi*i*f(k) 0]; + L = [0; 1]; + + % Solve A and Q + [A,Q] = lti_disc(F,L,Qc,dt); + + % The next state + x(:,k) = A*x(:,k-1) + chol(Q,'lower')*randn(size(Q,1),1); + + end + + z = z + H*x; + + end + + plot(t,5*(j-1)+z,'-k','LineWidth',.5) + + end + + xlabel('Time, $t$') + xlim([0 10.5]) + ylim([-3 18]) + + figure(2); clf + plot(t,f) + xlabel('Time, $t$') + ylabel('Frequency [Hz]') + + + \ No newline at end of file diff --git a/matlab/ch03_ex10_stochastic_spring_model.m b/matlab/ch03_ex10_stochastic_spring_model.m new file mode 100644 index 0000000..a1a9230 --- /dev/null +++ b/matlab/ch03_ex10_stochastic_spring_model.m @@ -0,0 +1,101 @@ +%% Example 3.10: Stochastic spring model +% +% Copyright: +% 2018 - Simo Särkkä and Arno Solin +% +% License: +% This software is provided under the MIT License. See the accompanying +% LICENSE file for details. + +%% +% Simulate model +% +if exist('rng') % Octave doesn't have rng + rng(1,'twister'); +else + randn('state',1); +end + +t1 = 10; +dt = 0.01; +steps = round(t1/dt); +g = 1; +v = 2; +q = 0.02; +x0 = [1;0]; +nmc = 10; + +F = [0 1; -v^2 -g]; +L = [0;1]; +[A,Q] = lti_disc(F,L,q,dt); + +M = []; +T = []; +m = x0; +t = 0; +for k=1:steps + m = A*m; + t = t + dt; + + M = [M m]; + T = [T t]; +end + +X1 = []; +X2 = []; +for n=1:nmc + X = []; + x = x0; + for k=1:steps + x = A*x + chol(Q,'lower')*randn(size(Q,1),1); + X = [X x]; + end + X1 = [X1;X(1,:)]; + X2 = [X2;X(2,:)]; +end + +figure(1); clf; hold on + plot(T,M(1,:),'-k','linewidth',1) + plot(T,X1','-','color',[.5 .5 .5]) + plot(T,M(1,:),'-k','linewidth',1) + legend('Deterministic solution','Realizations of the SDE'); + xlabel('Time, $t$'); ylabel('Displacement, $x(t)$') + box on + +%% +% Mean and covariance trajectories +% + +m = x0; +P = zeros(2,2); +MM = []; +VV = []; +for k=1:steps + m = A*m; + P = A*P*A' + Q; + + MM = [MM m]; + VV = [VV diag(P)]; +end + +figure(2); clf; hold on + + % Draw mean solution + plot(T,MM(1,:),'-k','linewidth',1); + + % Draw 95% shaded quantiles + ind = [1:10:steps steps]; + hp = fill([T(ind)'; flipud(T(ind)')], ... + [MM(1,ind)'-1.96*sqrt(VV(1,ind))'; flipud(MM(1,ind)'+1.96*sqrt(VV(1,ind))')],1); + set(hp,'EdgeColor',[.9 .9 .9],'FaceColor',[.9 .9 .9]) + + % Draw realization + plot(T,X1','-','color',[0.5 0.5 0.5]) + + % Draw mean on top of everything + plot(T,MM(1,:),'-k','linewidth',1); + + legend('Mean solution','95\% quantiles','Realizations of the SDE'); + xlabel('Time, $t$'); ylabel('Displacement, $x(t)$') + box on + set(gca,'Layer','Top') diff --git a/matlab/ch03_fig02_two_views_of_brownian_motion.m b/matlab/ch03_fig02_two_views_of_brownian_motion.m new file mode 100644 index 0000000..096627b --- /dev/null +++ b/matlab/ch03_fig02_two_views_of_brownian_motion.m @@ -0,0 +1,76 @@ +%% Figure 3.2: Two views of Brownian motion +% +% Copyright: +% 2018 - Simo Särkkä and Arno Solin +% +% License: +% This software is provided under the MIT License. See the accompanying +% LICENSE file for details. + +%% + + if exist('rng') % Octave doesn't have rng + rng(10,'twister') + else + randn('state',3); + end + + t = linspace(0,10,1000); + dt = t(2)-t(1); + x = [0 cumsum(sqrt(dt)*randn(1,numel(t)-1))]; + + % For Octave + norminv = @(x,m,s) -sqrt(2)*erfcinv(2*x) .* s + m; + normpdf = @(x,m,s) 1./(s .* sqrt(2*pi)) .* exp(-0.5*(x-m).^2 ./ s.^2); + + % Upper and lower bounds + ub = norminv(0.975,zeros(size(t)),sqrt(t)); + lb = norminv(0.025,zeros(size(t)),sqrt(t)); + + figure(1); clf + plot(t,x,'-k', ... + t,0*t,':k',... + t,ub,'--k', ... + t,lb,'-.k','LineWidth',0.5) + + % Legend + legend('Sample path of $\beta(t)$','Mean', ... + 'Upper 95\% quantile','Lower 95\% quantile') + + xlabel('Time, $t$') + ylabel('$\beta(t)$') + + +%% Mesh + + t = linspace(0,10,21); t(1) = .1; + x = linspace(-10,10,21); + [T,X] = meshgrid(t,x); + + P = normpdf(X(:),0.*T(:),sqrt(T(:))); + P = reshape(P,size(T)); + + figure(2); clf + + % Show surface + if exist('surf2patch') % Matlab + p=surf2patch(T,X,P); + h=patch(p,'EdgeColor','k','FaceColor',[.7 .7 .7],'FaceAlpha',.8,'LineWidth',.25); + + xlabel('Time, $t$') + ylabel('$\beta(t)$') + zlabel('$p(\beta(t))$') + + camproj('perspective') + else % Octave + surfl(T,X,P); + colormap gray; + end + + axis vis3d, grid on, box on + xlim([0 10]) + ylim([-10 10]) + zlim([0 1]) + + view([25 25]) + \ No newline at end of file diff --git a/matlab/ch03_fig10_white_noise.m b/matlab/ch03_fig10_white_noise.m new file mode 100644 index 0000000..adf58d4 --- /dev/null +++ b/matlab/ch03_fig10_white_noise.m @@ -0,0 +1,24 @@ +%% Figure 3.10: White noise +% +% Copyright: +% 2018 - Simo Särkkä and Arno Solin +% +% License: +% This software is provided under the MIT License. See the accompanying +% LICENSE file for details. + +%% A sketch of white noise + + if exist('rng') % Octave doesn't have rng + rng(0,'twister'); + else + randn('state',0); + end + + figure(1); clf + T = 0:0.001:1; + X = randn(size(T)); + plot(T,X); + xlabel('Time, $t$'); ylabel('$w(t)$') + + diff --git a/matlab/ch04_ex05_ou_process.m b/matlab/ch04_ex05_ou_process.m new file mode 100644 index 0000000..3d0e41f --- /dev/null +++ b/matlab/ch04_ex05_ou_process.m @@ -0,0 +1,58 @@ +%% Example 4.5: Solution of the Ornstein–Uhlenbeck process +% +% Copyright: +% 2018 - Simo Särkkä and Arno Solin +% +% License: +% This software is provided under the MIT License. See the accompanying +% LICENSE file for details. + +%% Simulate trajectories from the OU process + + if exist('rng') % Octave doesn't have rng + rng(10,'twister') + else + randn('state',1) + end + + lambda = 0.5; + q = 1; + dt = 0.01; + T = (0:dt:1); + x0 = 4; + M = exp(-lambda*T)*x0; + P = q/(2*lambda)*(1 - exp(-2*lambda*T)); + + XX = zeros(50,length(T)); + for n=1:size(XX,1) + x = x0; + for k=1:length(T) + XX(n,k) = x; + x = x - lambda * x * dt + sqrt(dt)*randn; + end + end + + + figure(1); clf; hold on + + % Shade the 95% quantiles + fill([T fliplr(T)],[M+1.96*sqrt(P) fliplr(M-1.96*sqrt(P))],1, ... + 'FaceColor',[.9 .9 .9],'EdgeColor',[.9 .9 .9]) + + % Plot realizations + h1 = plot(T,XX,'-','Color',[.5 .5 .5],'LineWidth',0.5); + + % Plot mean and quantiles + h2 = plot(T,M,'k-','LineWidth',1); + h34 = plot(T,M+1.96*sqrt(P),'--k', ... + T,M-1.96*sqrt(P),'--k','LineWidth',0.7); + + %set(h(1:3),'Linewidth',2) + + legend([h2(1) h34(1) h1(1)],'Mean','95\% quantiles','Realizations') + xlabel('Time, $t$'); ylabel('$x(t)$') + + ylim([0 5]) + box on + set(gca,'Layer','Top') + \ No newline at end of file diff --git a/matlab/ch04_fig01_brownian_motion.m b/matlab/ch04_fig01_brownian_motion.m new file mode 100644 index 0000000..fdc6446 --- /dev/null +++ b/matlab/ch04_fig01_brownian_motion.m @@ -0,0 +1,22 @@ +%% Figure 4.1: Brownian motion +% +% Copyright: +% 2018 - Simo Särkkä and Arno Solin +% +% License: +% This software is provided under the MIT License. See the accompanying +% LICENSE file for details. + +%% A sketch of Brownian motion + + if exist('rng') % Octave doesn't have rng + rng(4,'twister'); + else + randn('state',2); + end + + figure(1); clf + T = 0:0.001:1; + X = randn(size(T)); + plot([0 T],[0 cumsum(X)*sqrt((T(2)-T(1)))]); + xlabel('Time, $t$'); ylabel('$\beta(t)$') diff --git a/matlab/ch07_ex03_karhunen_loeve_series.m b/matlab/ch07_ex03_karhunen_loeve_series.m new file mode 100644 index 0000000..c863348 --- /dev/null +++ b/matlab/ch07_ex03_karhunen_loeve_series.m @@ -0,0 +1,48 @@ +%% Example 7.3: Karhunen–Loeve series of Brownian motion +% +% Copyright: +% 2018 - Simo Särkkä and Arno Solin +% +% License: +% This software is provided under the MIT License. See the accompanying +% LICENSE file for details. + +%% +% Series expansion +% + T = 10; + NN = [10 100 1000]; + + wd = [2.5 1.5 0.5]; + cl = [0.8 0.6 0.0]; + + legs = {}; + + figure(3); clf + for i=1:length(NN) + N = NN(i); + if exist('rng') % Octave doesn't have rng + rng(7,'twister') + else + randn('state',2) + end + n = 1:N; + t = 0:0.01:T; + lam = (2*T./(2*n-1)*pi).^2; + [tt,nn] = meshgrid(t,n); + phi = sqrt(2/T)*sin((2*nn-1)*pi.*tt./(2*T)); + z = sqrt(lam) .* randn(size(lam)); + h = plot(t,z*phi); + set(h,'LineWidth',wd(i)); + set(h,'Color',cl(i)*[1 1 1]); + hold on; + + legs{end+1} = sprintf('$N = %d$', NN(i)); + end + grid on; + + xlabel('Time, $t$') + ylabel('$\beta(t)$') + + legend(legs(:)); + \ No newline at end of file diff --git a/matlab/ch07_ex12_doobs_h_transform.m b/matlab/ch07_ex12_doobs_h_transform.m new file mode 100644 index 0000000..2c98515 --- /dev/null +++ b/matlab/ch07_ex12_doobs_h_transform.m @@ -0,0 +1,60 @@ +%% Example 7.12: Doob's h-transform +% +% Copyright: +% 2018 - Simo Särkkä and Arno Solin +% +% License: +% This software is provided under the MIT License. See the accompanying +% LICENSE file for details. + +%% For the OU model + + % Lock random seed + if exist('rng') % Octave doesn't have rng + rng(0,'twister'); + else + randn('state',2); + end + + % Parameters + lambda = 1; + q = 1; + x0 = 0; + T = 1; + xT = 5; + + % Define the derived functions + a = @(t) exp(-lambda*(T-t)); + sigma2 = @(t) q/(2*lambda) * (1-exp(-2*lambda*(T-t))); + + % Specify the model to simulate from + f = @(x,t) -lambda*x + q*a(t)/sigma2(t)*(xT - a(t)*x); + L = @(x,t) 1; + + % Simualte trajectories + t = linspace(0,T,100); + n = 100; + + % Allocate space + x = zeros(n,numel(t)); + + % Simulate using Euler-Maruyama + for i=1:n + x(i,:) = eulermaruyama(f,L,t,x0,q); + end + + figure(1); clf; hold on + + % Plot sample trajectories + plot(t,x','-','Color',[.5 .5 .5],'LineWidth',0.1) + + % Plot dashed line showing x(T) + plot([0 T],xT*[1 1],'--k') + + box on + xlabel('Time, $t$'); ylabel('$x(t)$') + set(gca,'Layer','Top') + set(gca,'YTick',[0:2:4 xT],'YTickLabel',{'$0$','$2$','$4$','$x_T$'}) + + + \ No newline at end of file diff --git a/matlab/ch07_ex17_feynman_kac.m b/matlab/ch07_ex17_feynman_kac.m new file mode 100644 index 0000000..1f3b361 --- /dev/null +++ b/matlab/ch07_ex17_feynman_kac.m @@ -0,0 +1,309 @@ +%% Example 7.17: Solution of an elliptic PDE using SDE simulation +% +% Copyright: +% 2018 - Simo Särkkä and Arno Solin +% +% License: +% This software is provided under the MIT License. See the accompanying +% LICENSE file for details. + +%% +% The equation to be considered is +% +% f_1(x,y) du/dx + f_2 du/dy +% + q/2 [d^2u(x,y)/dx^2 + d^2u(x,y)/dy^2] - r u(x,y) = 0 +% u(x,0) = psi_x1(x) +% u(x,L) = psi_x2(x) +% u(0,y) = psi_y1(y) +% u(L,y) = psi_y2(y) +% +% Finite differences approxmation +% +% The discretizations are +% +% d^2u(x,y)/dx^2 + d^2u(x,y)/dy^2 =~ (u(x+h,y) + u(x,y+h) - 4 u(x,y) + u(x-h,y) + u(x,y-h))/h^2 +% du/dx =~ (u(x+h,y) - u(x-h,y)) / (2h) +% du/dy =~ (u(x,y+h) - u(x,y-h)) / (2h) +% +% Giving +% +% f_1(x,y) (u(x+h,y) - u(x-h,y)) / (2h) +% + f_2(x,y) (u(x,y+h) - u(x,y-h)) / (2h) +% + q/2 (u(x+h,y) + u(x,y+h) - 4 u(x,y) + u(x-h,y) + u(x,y-h))/h^2 +% -r u(x,y) = 0 +% +% i.e. +% +% [q/2/h^2 + f_1(x,y)/(2h)] u(x+h,y) +% + [q/2/h^2 + f_2(x,y)/(2h)] u(x,y+h) +% - [2q/h^2 + r] u(x,y) +% + [q/2/h^2 - f_1(x,y)/(2h)] u(x-h,y) +% + [q/2/h^2 - f_2(x,y)/(2h)] u(x,y-h) = 0 +% +% with discretizations 0:h:(n+1)h the boundary conditions translate to +% +% when x=0: +% [q/2/h^2 + f_1(x,y)/(2h)] u(x+h,y) +% + [q/2/h^2 + f_2(x,y)/(2h)] u(x,y+h) +% - [2q/h^2 + r] u(x,y) +% + [q/2/h^2 - f_2(x,y)/(2h)] u(x,y-h) +% = -[q/2/h^2 - f_1(x,y)/(2h)] psi_y1(y) +% +% when x=L: +% [q/2/h^2 + f_2(x,y)/(2h)] u(x,y+h) +% - [2q/h^2 + r] u(x,y) +% + [q/2/h^2 - f_1(x,y)/(2h)] u(x-h,y) +% + [q/2/h^2 - f_2(x,y)/(2h)] u(x,y-h) +% = -[q/2/h^2 + f_1(x,y)/(2h)] psi_y2(y) +% +% when y=0: +% [q/2/h^2 + f_1(x,y)/(2h)] u(x+h,y) +% + [q/2/h^2 + f_2(x,y)/(2h)] u(x,y+h) +% - [2q/h^2 + r] u(x,y) +% + [q/2/h^2 - f_1(x,y)/(2h)] u(x-h,y) +% = -[q/2/h^2 - f_2(x,y)/(2h)] psi_x1(x) +% +% when y=L: +% [q/2/h^2 + f_1(x,y)/(2h)] u(x+h,y) +% - [2q/h^2 + r] u(x,y) +% + [q/2/h^2 - f_1(x,y)/(2h)] u(x-h,y) +% + [q/2/h^2 - f_2(x,y)/(2h)] u(x,y-h) +% = -[q/2/h^2 + f_2(x,y)/(2h)] psi_x2(x) +% +% which can be added to the system as additional equations. +% If we let +% +% U = [u(h,h) u(2h,h) ... u(nh,h) u(h,2h) u(2h,2h) ... u(nh,nh)]^T +% +% then we have for j,i=1,...,n +% +% u(ih,jh) = U((j-1)*n+i) +% + + % The model + lam1 = 0.1; + lam2 = 0.1; + L = 1; + f1 = @(x,y) -lam1*x; + f2 = @(x,y) -lam2*y; + psi_x1 = @(x) sin(2*pi*x); + psi_x2 = @(x) 1-x; + psi_y1 = @(y) y; + psi_y2 = @(y) sin(2*pi*y); + q = 0.1; + r = 1; + + % Discretization + n = 50; + h = L/(n+1); + + xgrid = h:h:n*h; + ygrid = h:h:n*h; + [YY,XX] = meshgrid(xgrid,ygrid); + X = XX(:); + Y = YY(:); + + % Evaluate boundary and f(x,y) + FF1 = f1(XX,YY); + FF2 = f2(XX,YY); + F1 = FF1(:); + F2 = FF2(:); + Px1 = psi_x1(xgrid); + Px2 = psi_x2(xgrid); + Py1 = psi_y1(ygrid); + Py2 = psi_y2(ygrid); + +% [q/2/h^2 + f_1(x,y)/(2h)] u(x+h,y) +% + [q/2/h^2 + f_2(x,y)/(2h)] u(x,y+h) +% - [2q/h^2 + r] u(x,y) +% + [q/2/h^2 - f_1(x,y)/(2h)] u(x-h,y) +% + [q/2/h^2 - f_2(x,y)/(2h)] u(x,y-h) + + F = spdiags((-2*q/h^2 - r)*ones(n*n,1),0,n*n,n*n); + b = spalloc(n*n,1,4*n); + for i=1:n + for j=1:n + if j~=1 + % u(x,y-h) + F((j-1)*n+i,(j-2)*n+i) = q/2/h^2 - F2((j-1)*n+i)/(2*h); + else + % boundary y=0 + b((j-1)*n+i) = b((j-1)*n+i) - (q/2/h^2 - F2((j-1)*n+i)/(2*h)) * Px1(i); + end + if j~=n + % u(x,y+h) + F((j-1)*n+i,(j-0)*n+i) = q/2/h^2 + F2((j-1)*n+i)/(2*h); + else + % boundary y=L + b((j-1)*n+i) = b((j-1)*n+i) - (q/2/h^2 + F2((j-1)*n+i)/(2*h)) * Px2(i); + end + if i~=1 + % u(x-h,y) + F((j-1)*n+i,(j-1)*n+i-1) = q/2/h^2 - F1((j-1)*n+i)/(2*h); + else + % boundary x=0 + b((j-1)*n+i) = b((j-1)*n+i) - (q/2/h^2 - F1((j-1)*n+i)/(2*h)) * Py1(j); + end + if i~=n + % u(x+h,y) + F((j-1)*n+i,(j-1)*n+i+1) = q/2/h^2 + F1((j-1)*n+i)/(2*h); + else + % boundary x=L + b((j-1)*n+i) = b((j-1)*n+i) - (q/2/h^2 + F1((j-1)*n+i)/(2*h)) * Py2(j); + end + end + end + +% clf; +% subplot(2,1,1); +% spy(F); + + U = F\b; + UU = reshape(full(U),size(XX)); + + figure(1); clf; + surf(XX,YY,UU); + shading interp; + if exist('camproj') % Octave doesn't have this + camproj perspective; + end + colormap winter; + colorbar; + + xlabel('x'); + ylabel('y'); + +%% +% Feynman-Kac approximation +% + % Lock random seed + if exist('rng') % Octave doesn't have rng + rng(1); + else + randn('state',2); + end + + dt = 0.01; + a1 = exp(-lam1*dt); + a2 = exp(-lam2*dt); + sq1 = sqrt(q/(2*lam1)*(1 - exp(-2*lam1*dt))); + sq2 = sqrt(q/(2*lam2)*(1 - exp(-2*lam2*dt))); + + UFK = zeros(size(XX)); + NN = zeros(size(XX)); + + squ_x = [0 L L 0]; + squ_y = [0 0 L L]; + + stored_SX1 = []; + stored_SX2 = []; + stored_whos_in = []; + + nmc = 100; + for mc=1:nmc + SX1 = XX(:); + SX2 = YY(:); + + count = 0; + whos_in = true(1,length(SX1)); + num_in = sum(whos_in); + while num_in > 0 + SX1(whos_in) = a1 * SX1(whos_in) + sq1 * randn(size(SX1(whos_in))); + SX2(whos_in) = a2 * SX2(whos_in) + sq2 * randn(size(SX2(whos_in))); + + ind = find(whos_in(:) & (SX1(:) < 0)); + whos_in(ind) = false; + UFK(ind) = UFK(ind) + psi_y1(SX2(ind)); + NN(ind) = NN(ind) + 1; + + ind = find(whos_in(:) & (SX1(:) > L)); + whos_in(ind) = false; + UFK(ind) = UFK(ind) + psi_y2(SX2(ind)); + NN(ind) = NN(ind) + 1; + + ind = find(whos_in(:) & (SX2(:) < 0)); + whos_in(ind) = false; + UFK(ind) = UFK(ind) + psi_x1(SX1(ind)); + NN(ind) = NN(ind) + 1; + + ind = find(whos_in(:) & (SX2(:) > L)); + whos_in(ind) = false; + UFK(ind) = UFK(ind) + psi_x2(SX1(ind)); + NN(ind) = NN(ind) + 1; + + num_in = sum(whos_in); + + count = count + 1; + + if mc == 1 && count == 200 + stored_whos_in = whos_in; + stored_SX1 = SX1; + stored_SX2 = SX2; + clf; + plot(SX1(:),SX2(:),'.',squ_x,squ_y,'k'); + pause(1); + end + + if rem(count,100) == 0 + + subplot(2,2,1); + pcolor(XX,YY,UU); + colorbar; + + subplot(2,2,2); + pcolor(XX,YY,UFK ./ NN); + colorbar; + + subplot(2,2,3); + plot(SX1(:),SX2(:),'.',squ_x,squ_y,'k'); + + subplot(2,2,4); + pcolor(XX,YY,abs(UU - UFK ./ NN)); + title(sprintf('MC iteration %d/%d',mc,nmc)); + colorbar; + + drawnow; + fprintf('count = %d, num_in = %d\n',count,num_in); + end + + end + end + + % Consider caching the results for later use + % save fkac1 + +%% +% Plot the results +% + figure(1); clf; + imagesc(XX(:,1),YY(1,:),UU); + caxis([-1 1]) + shading flat + colormap gray + colorbar + axis xy + + figure(2); clf; + imagesc(XX(:,1),YY(1,:),UFK ./ NN); + caxis([-1 1]) + shading flat + colormap gray + colorbar; + axis xy + + figure(3); clf; + plot(stored_SX1(stored_whos_in),stored_SX2(stored_whos_in),'kx',... + stored_SX1(~stored_whos_in),stored_SX2(~stored_whos_in),'ko',squ_x,squ_y,'k'); + axis xy + + figure(4); clf; + imagesc(XX(:,1),YY(1,:),abs(UU - UFK ./ NN)); + shading flat + colormap gray + colorbar + axis xy + + nmc + + err = UU - UFK ./ NN; + rmse = mean(err(:).^2) \ No newline at end of file diff --git a/matlab/ch08_ex06_weak_itotaylor.m b/matlab/ch08_ex06_weak_itotaylor.m new file mode 100644 index 0000000..5073fc7 --- /dev/null +++ b/matlab/ch08_ex06_weak_itotaylor.m @@ -0,0 +1,169 @@ +%% Example 8.6: Simulating from a trigonometric nonlinear SDE +% +% Copyright: +% 2018 - Simo Särkkä and Arno Solin +% +% License: +% This software is provided under the MIT License. See the accompanying +% LICENSE file for details. + +%% Gaussian approximation + + % Lock random seed + if exist('rng') % Octave doesn't have rng + rng(1,'twister') + else + randn('state',1); + rand('state',1); + end + + % Parameters + tspan = 0:1:10; + + % The model + f = @(x,t) -(1/10)^2*sin(x).*cos(x).^3; + L = @(x,t) 1/10*cos(x).^2; + + % The derivatives + df = @(x,t) -1/100*cos(x).^2.*(2*cos(2*x)-1); + ddf = @(x,t) 1/100*(sin(2*x) + 2*sin(4*x)); + dL = @(x,t) -1/5*sin(x).*cos(x); + ddL = @(x,t) -1/5*cos(2*x); + + % Exact solution + solfun = @(wt,x0) atan(1/10*wt +tan(x0)); + + % Intial + x0 = 1; + + +%% Sample + + x = zeros(1,10000); + xem = x; + xw20 = x; + xw20g = x; + + % Lock random seed + if exist('rng') % Octave doesn't have rng + rng(1,'twister') + else + randn('state',1); + rand('state',1); + end + + for j=1:size(x,2) + + % Use Euler-Maruyama + foo = eulermaruyama_weak(f,L,tspan,x0,1); + xem(:,j) = foo(:,end); + + % Report + if rem(j,100)==0, j, end + + end + + % Lock random seed + if exist('rng') % Octave doesn't have rng + rng(1,'twister') + else + randn('state',1); + rand('state',1); + end + + for j=1:size(x,2) + + % Weak order 2.0 + foo = w20scalar({f,df,ddf},{L,dL,ddL},tspan,x0,1,false); + xw20(:,j) = foo(:,end); + + % Report + if rem(j,100)==0, j, end + + end + + % Lock random seed + if exist('rng') % Octave doesn't have rng + rng(1,'twister') + else + randn('state',1); + rand('state',1); + end + + for j=1:size(x,2) + + % Weak order 2.0 (Gaussian increments) + foo = w20scalar({f,df,ddf},{L,dL,ddL},tspan,x0,1,true); + xw20g(:,j) = foo(:,end); + + % Store + x(:,j) = foo(:,end); + + % Report + if rem(j,100)==0, j, end + + end + + +%% Visualize + + % Lock random seed + if exist('rng') % Octave doesn't have rng + rng(1,'twister') + else + randn('state',1); + rand('state',1); + end + + % Samples from the exact solution + xe = solfun(sqrt(tspan(end))*randn(1,200000),x0); + + % Bins + nbins = 64; + t = linspace(min(xe),max(xe)+.1,nbins); + + figure(2); clf; hold on + + % Show solution + n = histc(xe,t); + fill(t,n/numel(xe),1,'FaceColor',[.7 .7 .7],'EdgeColor',[.7 .7 .7]) + + % Show solution w2.0 + n = histc(xw20,t); + stairs(t-(t(2)-t(1))/2,n/numel(xw20),'-k') + + % Limits + xlim([.35 1.25]) + lims = ylim; + + % Label + xlabel('$x$') + + % Ticks + set(gca,'XTick',0:.2:1.2) + + + figure(3); clf; hold on + + % Show solution + n = histc(xe,t); + fill(t,n/numel(xe),1,'FaceColor',[.7 .7 .7],'EdgeColor',[.7 .7 .7]) + + % Show solution w2.0 (Gaussian increments) + n = histc(xw20g,t); + stairs(t-(t(2)-t(1))/2,n/numel(xw20g),'-k') + + % Limits + %xlim([min(t) max(t)]) + xlim([.35 1.25]) + ylim(lims) + + % Show legend + legend('Exact','Weak order $2.0$') + + % Label + xlabel('$x$') + + % Ticks + set(gca,'XTick',0:.2:1.2) + diff --git a/matlab/ch08_ex11_ode_solvers.m b/matlab/ch08_ex11_ode_solvers.m new file mode 100644 index 0000000..d92f0e5 --- /dev/null +++ b/matlab/ch08_ex11_ode_solvers.m @@ -0,0 +1,113 @@ +%% Example 8.11: Comparison of ODE solvers +% +% Copyright: +% 2018 - Simo Särkkä and Arno Solin +% +% License: +% This software is provided under the MIT License. See the accompanying +% LICENSE file for details. + +%% Plot phase portrait + + % Define the system in terms of the differential equation + f = @(x) [x(1,:)-x(2,:)-x(1,:).^3; + x(1,:)+x(2,:)-x(2,:).^3]; + + % Define arrow + arrow1 = [-1 1 0 -1; -.5 -.5 2 -.5]'; + + % Define arrow without head + arrow2 = .4*[-1 -1 1 1 -1;-.1 .1 .1 -.1 -.1]'; + + % Phase portrait figure + figure(1); clf; hold on + + % Discrete grid + x = linspace(-2,2,15); + y = linspace(-2,2,15); + [X Y] = meshgrid(x,y); + + % Form vector field + UV = f([X(:) Y(:)]'); + U = reshape(UV(1,:),size(X)); + V = reshape(UV(2,:),size(Y)); + + % Plot field + X(X==0 & Y==0) = nan; + newquiver(X,Y,U,V,'X',arrow2, ... + 'color',[.7 .7 .7],'EdgeColor',[.7 .7 .7]) + + % Change axis + axis equal; axis([-2.2 2.2 -2.2 2.2]); + set(gca,'XTick',[-2:4:2],'YTick',[-2:4:2]) + + % Show some trajectories (start from outside) + x0 = 3*[cos(linspace(0,2*pi,25)); sin(linspace(0,2*pi,25))]; + + tin = 0:2^-4:10; + for i=1:6 + try + xout = euler(@(x,t) f(x),tin,x0(:,i)); xout = xout'; + plot(xout(:,1),xout(:,2),'-k','LineWidth',.5,'Color',[.5 .5 .5]) + uv = f(xout'); ind = 4; + newquiver(xout(ind,1),xout(ind,2),uv(2,ind),-uv(1,ind), ... + 'X',arrow1,'scale',.04) + end + end + for i=7:12 + try + xout = heun(@(x,t) f(x),tin,x0(:,i)); xout = xout'; + plot(xout(:,1),xout(:,2),'-k','LineWidth',.5) + uv = f(xout'); ind = 4; + newquiver(xout(ind,1),xout(ind,2),uv(2,ind),-uv(1,ind), ... + 'X',arrow1,'scale',.04) + end + end + for i=13:18 + try + xout = impliciteuler(@(x,t) f(x),tin,x0(:,i)); xout = xout'; + plot(xout(:,1),xout(:,2),'-k','LineWidth',.5,'Color',[.5 .5 .5]) + uv = f(xout'); ind = 4; + newquiver(xout(ind,1),xout(ind,2),uv(2,ind),-uv(1,ind), ... + 'X',arrow1,'scale',.04) + end + end + for i=19:24 + try + xout = rk4simple(@(x,t) f(x),tin,x0(:,i)); xout = xout'; + plot(xout(:,1),xout(:,2),'-k','LineWidth',.5) + uv = f(xout'); ind = 4; + newquiver(xout(ind,1),xout(ind,2),uv(2,ind),-uv(1,ind), ... + 'X',arrow1,'scale',.04) + end + end + + % Labels + text( 1.9, 1.9,'Forward Euler','HorizontalAlignment','right') + text(-1.9, 1.9,'Heun','HorizontalAlignment','left') + text(-1.9,-1.9,'Backward Euler','HorizontalAlignment','left') + text( 1.9,-1.9,'RK4','HorizontalAlignment','right') + + % Show some trajectories (start from inside) + x0 = .02*[cos(linspace(0,2*pi,5)); sin(linspace(0,2*pi,5))]; + for i=1:size(x0,2) + try + xout = rk4simple(@(x,t) f(x),linspace(0,10,48),x0(:,i)); xout = xout'; + plot(xout(:,1),xout(:,2),'-k','LineWidth',.5) + uv = f(xout'); ind = 18; + newquiver(xout(ind,1),xout(ind,2),uv(2,ind),-uv(1,ind), ... + 'X',arrow1,'scale',.04) + end + end + + % Axis labels + xlabel('$x_1$') + ylabel('$x_2$') + + % Plot fixed point + plot(0,0,'ok','MarkerFaceColor','w','MarkerSize',3) + + % Set figure options + set(gcf,'Color','w') + + \ No newline at end of file diff --git a/matlab/ch08_ex15_duffing_van_der_pol.m b/matlab/ch08_ex15_duffing_van_der_pol.m new file mode 100644 index 0000000..3f80640 --- /dev/null +++ b/matlab/ch08_ex15_duffing_van_der_pol.m @@ -0,0 +1,170 @@ +%% Example 8.15: Duffing van der Pol oscillator +% +% Copyright: +% 2018 - Simo Särkkä and Arno Solin +% +% License: +% This software is provided under the MIT License. See the accompanying +% LICENSE file for details. + +%% Duffing van der Pol + + % Time-span + tspan = 0:2^-5:20; + + % Parameters + alpha = 1; + + % Define arrow (for visalization) + arrow1 = [-1 1 0 -1; -.5 -.5 2 -.5]'; + + % The model + f = @(x,t) [x(2,:); x(1,:).*(alpha - x(1,:).^2)-x(2,:)]; + L = @(x,t) [zeros(1,size(x,2)); x(1,:)]; + + +%% ODE + + figure(1); clf; hold on + + for j=1:10 + %x = rk4(f,tspan,[-2-.2*j; 0]); + x = rk4simple(f,tspan,[-2-.2*j; 0]); + %[~,x] = ode45(@(t,x) f(x,t),tspan,[-2-.2*j; 0]); x = x'; + + % Plot trajectory + plot(x(1,:),x(2,:),'-k','LineWidth',.25) + + % Plot direction + uv = f(x,[]); ind = 10; + newquiver(x(1,ind),x(2,ind),uv(2,ind),-uv(1,ind), ... + 'X',arrow1,'scale',.04*[1 14.8/8.8]) + + end + + % Set axis limits + axis([-4.4 4.4 -4.8 10]), axis square + %set(gca,'XTick',-4:2:4,'YTick',-4:2:8) + xlabel('$x_1$'); ylabel('$x_2$') + box on + +%% SDE + + figure(2); clf; hold on + figure(3); clf; hold on + + for j=1:10 + + figure(2); + + % Lock random seed + if exist('rng') % Octave doesn't have rng + rng(3,'twister') + else + randn('state',2); + rand('state',2); + end + + % Use the strong order 1.0 method + x = srkS10scalarnoise(f,L,tspan,[-2-.2*j; 0],.5^2); + + % Plot trajectory + plot(x(1,:),x(2,:),'-k','LineWidth',.25) + + % Plot direction + uv = f(x,[]); ind = 10; + newquiver(x(1,ind),x(2,ind),uv(2,ind),-uv(1,ind), ... + 'X',arrow1,'scale',.04*[1 14.8/8.8]) + + figure(3); + plot(tspan,x(1,:),'-k') + plot(tspan,x(2,:),'-','Color',[.7 .7 .7]) + + end + + % Set axis limits + figure(2) + axis([-4.4 4.4 -4.8 10]), axis square + %set(gca,'XTick',-4:2:4,'YTick',-4:2:8) + xlabel('$x_1$'); ylabel('$x_2$') + box on + + % Set axis limits + figure(3) + xlim([0 20]) + xlabel('Time, $t$'); ylabel('$x$') + legend('$x_1(t)$','$x_2(t)$') + box on + + +%% Weak approximation + + % Time-span + tspan = 0:2^-4:20; + + % Reset random seed + % Lock random seed + if exist('rng') % Octave doesn't have rng + rng(1,'twister') + else + randn('state',2); + rand('state',2); + end + + % Initial point + x0 = [-3;0]; + + % Number of smaples + x = zeros(2,10000); + + % Iterate + for j=1:size(x,2) + + % Weak SRK scheme + foo = srkW20(f,L,tspan,x0,.5^2,true); + + % Store + x(:,j) = foo(:,end); + + % Report + if rem(j,100)==0, + figure(4); clf + hist(x(1,1:j),ceil(sqrt(j))) + drawnow + j + end + + end + + +%% Histogram + + % Bins for histogram + t = linspace(min(x(1,:)),max(x(1,:)),64); + + figure(5); clf + + % Show solution w2.0 + n = histc(x(1,:),t); + stairs(t-(t(2)-t(1))/2,n/size(x,2),'-k') + + % Label + xlabel('$x_1$') + + % Ticks + box off + xlim([-2.2 2.2]) + %set(gca,'XTick',0:.2:1.2) + + figure(6); clf + + % Show solution w2.0 + plot(x(1,:),x(2,:),'.k') + + % Label + xlabel('$x_1$') + ylabel('$x_2$') + + % Ticks + box on + \ No newline at end of file diff --git a/matlab/ch08_ex21_leapfrog_verlet.m b/matlab/ch08_ex21_leapfrog_verlet.m new file mode 100644 index 0000000..ef726cb --- /dev/null +++ b/matlab/ch08_ex21_leapfrog_verlet.m @@ -0,0 +1,181 @@ +%% Example +% +% Copyright: +% 2018 - Simo Särkkä and Arno Solin +% +% License: +% This software is provided under the MIT License. See the accompanying +% LICENSE file for details. + +%% Leapfrog / Verlet integration + +dts = logspace(-2,0,10); +err = nan(10,numel(dts)); + +for j=1:numel(dts) + fprintf('Running step %d/%d\n',j,numel(dts)); + +% Leapfrog / Verlet integration + + % Lock seed + if exist('rng') % Octave doesn't have rng + rng(j,'twister') + else + randn('state',j); + end + + % Parameters + n = 25000; + g = 1; % g>0 + eta = 1; + q = 1; + x0 = 1; + v0 = 0; + + % Time discretization + dt = dts(j); %0.1; + t = 0:dt:10; + + % Specify the model + f = @(x) -g*x; + s = @(x) 1; + + % Allocate space and set initial conditions + x = zeros(n,numel(t)); x(:,1) = x0; + v = zeros(n,numel(t)); v(:,1) = v0; + + % For each trajectory + for i=1:n + + % Simulate one trajectory from the leapfrog method + z = leapfrog(f,s,eta,t,[x0; v0],q); + + x(i,:) = z(1,:); + v(i,:) = z(2,:); + + end + + figure(1); clf + subplot(211) + plot(t,x(1:5,:), ... + t,mean(x), ... + t,mean(x)+1.96*std(x),'--', ... + t,mean(x)-1.96*std(x),'--') + xlabel('Time, t'), ylabel('Position, x(t)') + subplot(212) + plot(t,v(1:5,:), ... + t,mean(v), ... + t,mean(v)+1.96*std(v),'--', ... + t,mean(v)-1.96*std(v),'--') + xlabel('Time, t'), ylabel('Velocity, v(t)') + + +% Euler-Maruyama + + % Same random seed + % Lock seed + if exist('rng') % Octave doesn't have rng + rng(j,'twister') + else + randn('state',j); + end + + % Specify the model + f = @(x,t) [0 1; -g -eta]*x; + L = @(x,t) [0; 1]; + + % Allocate space and set initial conditions + x_em = zeros(n,numel(t)); x_em(:,1) = x0; + v_em = zeros(n,numel(t)); v_em(:,1) = v0; + + % For each trajectory + for i=1:n + + % Simulate one trajectory from the leapfrog method + z = eulermaruyama(f,L,t,[x0; v0],q); + + x_em(i,:) = z(1,:); + v_em(i,:) = z(2,:); + + end + + figure(1); clf + subplot(211) + plot(t,x_em(1:5,:), ... + t,mean(x_em), ... + t,mean(x_em)+1.96*std(x_em),'--', ... + t,mean(x_em)-1.96*std(x_em),'--') + xlabel('Time, t'), ylabel('Position, x(t)') + subplot(212) + plot(t,v(1:5,:), ... + t,mean(v_em), ... + t,mean(v_em)+1.96*std(v_em),'--', ... + t,mean(v_em)-1.96*std(v_em),'--') + xlabel('Time, t'), ylabel('Velocity, v(t)') + + +% Closed-form solution + + % This is a LTI SDE of the form + F = [0 1; -g -eta]; + L = [0; 1]; + Qc = q; + + % Discretize + [A,Q] = lti_disc(F,L,Qc,dt); + + % Evaluate + M = zeros(2,numel(t)); M(:,1) = [x0; v0]; + P = zeros(2,2,numel(t)); + + % Loop through + for k=1:numel(t)-1 + M(:,k+1) = A*M(:,k); + P(:,:,k+1) = A*P(:,:,k)*A'+Q; + end + + figure(2); clf + subplot(211) + plot(t,M(1,:), ... + t,M(1,:)'+1.96*sqrt(squeeze(P(1,1,:))),'--', ... + t,M(1,:)'-1.96*sqrt(squeeze(P(1,1,:))),'--') + subplot(212) + plot(t,M(2,:), ... + t,M(2,:)'+1.96*sqrt(squeeze(P(2,2,:))),'--', ... + t,M(2,:)'-1.96*sqrt(squeeze(P(2,2,:))),'--') + + +% Calculate errors + + err(1,j) = P(1,1,end)-std(x(:,end))^2; + err(2,j) = P(2,2,end)-std(v(:,end))^2; + + err(3,j) = P(1,1,end)-std(x_em(:,end))^2; + err(4,j) = P(2,2,end)-std(v_em(:,end))^2; + + err(5,j) = M(1,end)-mean(x(:,end)); + err(6,j) = M(2,end)-mean(v(:,end)); + + err(7,j) = M(1,end)-mean(x_em(:,end)); + err(8,j) = M(2,end)-mean(v_em(:,end)); + + err(9,j) = mean(abs(squeeze(P(1,1,:))'-std(x).^2)); + err(10,j) = mean(abs(squeeze(P(1,1,:))'-std(x_em).^2)); + + +end + + +%% Compose results + + figure(1); clf + loglog(dts,abs(err(9,:)),'-k', ... + dts,abs(err(10,:)),'--k') + xlabel('Time step length, $\Delta t$') + ylabel('Mean absolute error in $\sigma_x^2$') + legend('Leapfrog Verlet', 'Euler--Maruyama') + %ylim([0.4 10]) + xlim([min(dts) max(dts)]) + set(gca,'XTick',[0.01 0.1 1],'XTickLabel',[0.01 0.1 1]) + set(gca,'YTick',[1e-3 1e-2 1e-1 1e0 1e1], ... + 'YTickLabel',{'','$10^{-2}$','$10^{-1}$','$10^{0}$','$10^1$'}) diff --git a/matlab/ch08_ex24_exact_simulation.m b/matlab/ch08_ex24_exact_simulation.m new file mode 100644 index 0000000..fcffed3 --- /dev/null +++ b/matlab/ch08_ex24_exact_simulation.m @@ -0,0 +1,134 @@ +%% Example 8.24: Exact simulation of sine diffusion (EA1) +% +% Copyright: +% 2018 - Simo Särkkä and Arno Solin +% +% License: +% This software is provided under the MIT License. See the accompanying +% LICENSE file for details. + +%% + + % Parameters + T = 1.5; + n = 100000; + + % Define the SDE model + f = @(x,t) sin(x); + L = @(x,t) 1; + + % Define d/dx f(x) + df = @(x) cos(x); + + % Define psi(x) = int_0^x sin(u) du + psi = @(x) 1-cos(x); + + % Bounds: k1 < .5*(f(x)^2 + df(x)) < k2 + k1 = -0.5; + k2 = 0.6250; + M = (k2-k1); + + % Initialize h(x) = exp(psi(x)?x2/2T)/c + c = integral(@(u) exp(psi(u)-u.^2/(2*T)),-inf,inf); + h = @(u) 1/c*exp(psi(u)-u.^2/(2*T)); + + % Initialize: phi(x) = .5*f(x) + .5*df(x)?k1 + phi = @(x) 1/2*(f(x,0).^2+df(x))-k1; + + % Choose an epsilon for rejection sampling + epsilon = 1/70; + +%% + + % Lock seed + if exist('rng') % Octave doesn't have rng + rng(0,'twister') + else + randn('state',0); + rand('state',0); + end + + +% Allocate space for samples +xea = zeros(n,1); + +for i=1:n + + % Proposal skeleton + while (true) + + % Step #1 + lambda = T*M; + + if exist('poissrnd') + tau = poissrnd(lambda); + else + % Knuth's algorithm for generating tau = Poisson(lambda) + ell = exp(-lambda); + tau = 0; + p = 1; + while p > ell + tau = tau + 1; + u = rand; + p = p * u; + end + tau = tau - 1; + end + + % Step #2 + t = T*rand(1,tau); + x = M*rand(1,tau); + + % Step #3 + beta = zeros(1,2+tau); + beta(1) = 0; + beta(2) = rejectionsample(h,epsilon); + tb = [0 T t]; + + % The rest of the points by Brownian bridges + for j=3:numel(beta) + + % Lower + ind = find(tb(1:j-1)tb(j)); + [t2,i2] = min(tb(ind)); + b2 = beta(ind); b2 = b2(i2); + + % The Brownian bridge + mu = b1+(b2-b1)/(t2-t1)*(tb(j)-t1); + sigma2 = (t2-tb(j))*(tb(j)-t1)/(t2-t1); + beta(j) = mu + sqrt(sigma2)*randn(1); + + end + + % Step #4 + N = sum(x= t_k) + % + % p(x | Y_k) = + % 1/sqrt{2 pi P} exp(-(x-m)^2/(2P)) exp(-1/2 P) cosh(x) / cosh(m) + % + % and + % + % E[x] = m + tanh(m) P + % E[x^2] = m^2 + P^2 + 2 m tanh(m) P + P + % V[x] = P + [1 - tanh(m)^2] P^2 + + dx = 0.001; + xx = -15:dx:15; + dens = @(x,m,P) 1/sqrt(2*pi*P)*exp(-P/2)/cosh(m) * exp(-(x-m).^2/(2*P)) .* cosh(x); + + m = 1; + P = 5; + pp = dens(xx,m,P); + figure(1); clf; plot(xx,pp); + + sum(pp .* dx) + [sum(xx .* pp * dx) m+tanh(m)*P] + [sum(xx.^2 .* pp * dx) m^2+P^2+2*m*tanh(m)*P+P] + [sum((xx-m-tanh(m)*P).^2 .* pp * dx) P+(1-tanh(m)^2)*P^2] + +%% +% Generate data +% + % Lock random seed + if exist('rng') % Octave doesn't have rng + rng(1) + else + randn('state',9); + end + + steps = 500; + meas_mod = 20; + dt = 0.01; + R = 0.1; + + X = []; + Y = []; + T = []; + x0 = 0; + x = x0; + t = 0; + for k=1:steps + x = x + tanh(x) * dt + sqrt(dt)*randn; + if rem(k,meas_mod) == 0 + y = x + sqrt(R)*randn; + else + y = NaN; + end + t = t + dt; + X(k) = x; + Y(k) = y; + T(k) = t; + end + + figure(1); clf + plot(T,X,T,Y,'o'); + +%% +% Benes-Daum filter +% + EB = zeros(1,length(Y)); + VB = zeros(1,length(Y)); + MMB = zeros(1,length(Y)); + PPB = zeros(1,length(Y)); + + mb0 = 0.1; + Pb0 = 1; + + mb = mb0; + Pb = Pb0; + + m0 = mb + Pb * tanh(mb); + P0 = Pb + (1-tanh(mb)^2) * Pb^2; + + for k=1:length(Y) + % + % The Benes-Daum filter + % + Pb = Pb + dt; + if ~isnan(Y(k)) + mb = mb + Pb / (R + Pb)*(Y(k) - mb); + Pb = Pb - Pb^2/(R + Pb); + end + MMB(:,k) = mb; + PPB(:,k) = Pb; + EB(:,k) = mb + Pb * tanh(mb); + VB(:,k) = Pb + (1-tanh(mb)^2) * Pb^2; + end + + figure(1); clf; + subplot(2,1,1); + plot(T,MMB,T,EB); + subplot(2,1,2); + plot(T,PPB,T,VB); + + rmse_b = sqrt(mean((EB - X).^2)) + +%% +% Determine the quantiles of Benes solution +% + QQ1 = zeros(size(MMB)); + QQ2 = zeros(size(MMB)); + for k=1:length(MMB) + pp = dens(xx,MMB(:,k),PPB(:,k)); + cu = cumsum(pp ./ sum(pp)); + ind1 = find(cu >= 2.5/100,1); + ind2 = find(cu >= 97.5/100,1); + QQ1(k) = xx(ind1); + QQ2(k) = xx(ind2); + end + + figure(1); clf; + plot(T,EB,T,QQ1,'--',T,QQ2,'--'); + +%% +% EKF solution to the same problem +% + me = m0; + Pe = P0; + MMe = zeros(1,length(Y)); + PPe = zeros(1,length(Y)); + ekf_steps = 5; + + for k=1:length(Y) + dt2 = dt/ekf_steps; + for i=1:ekf_steps + Pe = Pe + dt2 * (2 * (1 - tanh(me)^2) * Pe + 1); + me = me + dt2 * tanh(me); + end + if ~isnan(Y(k)) + me = me + Pe / (R + Pe) * (Y(k) - me); + Pe = Pe - Pe^2/(R + Pe); + end + MME(:,k) = me; + PPE(:,k) = Pe; + end + + rmse_e = sqrt(mean((MME - X).^2)) + + figure(1); clf; + subplot(2,2,1); + plot(T,EB,T,QQ1,'--',T,QQ2,'--'); + subplot(2,2,2); + plot(T,MME,T,MME-1.96*sqrt(PPE),'--',T,MME+1.96*sqrt(PPE),'--'); + + subplot(2,2,3); + plot(T,EB,T,MME) + + subplot(2,2,4); + plot(T,VB,T,PPE) + +%% +% ERTS Smoother +% + MMS = MME; + PPS = PPE; + ms = MME(end); + Ps = PPE(end); + for k=length(Y)-1:-1:1 + m = MME(k); + P = PPE(k); + dms = dt * tanh(m) + dt * ((1 - tanh(m)^2) * P + 1)/P * (ms - m); + dPs = dt * 2 * ((1 - tanh(m)^2) * P + 1)/P * Ps - dt; + ms = ms - dms; + Ps = Ps - dPs; + MMS(k) = ms; + PPS(k) = Ps; + end + + rmse_s = sqrt(mean((MMS - X).^2)) + + figure(1); clf; + plot(T,MMS,T,MMS-1.96*sqrt(PPS),'--',T,MMS+1.96*sqrt(PPS),'--'); + +%% +% Plot the final figures +% + + % + % Benes-Daum + % + figure(1); clf; hold on + + % Plot mean + h0 = plot(T,EB,'-k'); + set(h0,'LineWidth',1); + set(h0,'Color',[0.5 0.5 0.5]); + + % Plot uncetainty + h1 = fill([T'; flipud(T')], [QQ1'; flipud(QQ2')],1); + set(h1,'EdgeColor',[.7 .7 .7],'FaceColor',[.7 .7 .7]) + + % Plot signal + h2 = plot(T,X,'-k'); + set(h2,'LineWidth',1); + + % Plot mean + h3 = plot(T,EB,'-k'); + set(h3,'LineWidth',1); + set(h3,'Color',[0.5 0.5 0.5]); + + % Plot observations + h4 = plot(T,Y,'+k'); + + box on + xlabel('Time, $t$'), ylabel('Signal, $x(t)$') + + legend([h0, h1, h2, h4],'BD mean','BD 95\% quantiles','Signal','Observations') + + % + % Benes-Daum vs EKF + % + figure(2); clf; hold on + + % Plot BD mean + h1 = plot(T,EB,'-k'); + set(h1,'LineWidth',1); + set(h1,'Color',[0.5 0.5 0.5]); + + % Plot EKF mean + h2 = plot(T,MME,'--k'); + set(h2,'LineWidth',1); + %set(h2,'Color',[0.0 0.0 0.0]); + + % Plot BD uncertainty + h3 = fill([T'; flipud(T')], [QQ1'; flipud(QQ2')],1); + set(h3,'EdgeColor',[.7 .7 .7],'FaceColor',[.7 .7 .7]) + + % Plot BD mean + h4 = plot(T,EB,'-k'); + set(h4,'LineWidth',1); + set(h4,'Color',[0.5 0.5 0.5]); + + % Plot EKF mean + h5 = plot(T,MME,'--k'); + set(h5,'LineWidth',1); + set(h5,'Color',[0.0 0.0 0.0]); + + % Plot EKF uncertainty + h6 = plot(T,MME-1.96*sqrt(PPE),'-.k',T,MME+1.96*sqrt(PPE),'-.k'); + set(h6,'LineWidth',.5); + + box on + xlabel('Time, $t$'), ylabel('Signal, $x(t)$') + legend([h1, h2, h3, h6(1)],'BD mean','EKF mean','BD 95\% quantiles','EKF 95\% quantiles') + + % + % ERTS + % + figure(3); clf; hold on + + % Plot mean + h0 = plot(T,MMS,'-k'); + set(h0,'LineWidth',1); + set(h0,'Color',[0.5 0.5 0.5]); + + QQ1S = MMS-1.96*sqrt(PPS); + QQ2S = MMS+1.96*sqrt(PPS); + + % Plot uncetainty + h1 = fill([T'; flipud(T')], [QQ1S'; flipud(QQ2S')],1); + set(h1,'EdgeColor',[.7 .7 .7],'FaceColor',[.7 .7 .7]) + + % Plot signal + h2 = plot(T,X,'-k'); + set(h2,'LineWidth',1); + + % Plot mean + h3 = plot(T,MMS,'-k'); + set(h3,'LineWidth',1); + set(h3,'Color',[0.5 0.5 0.5]); + + % Plot observations + h4 = plot(T,Y,'+k'); + + box on + xlabel('Time, $t$'), ylabel('Signal, $x(t)$') + + legend([h0, h1, h2, h4],'ERTS mean','ERTS 95\% quantiles','Signal','Observations') + diff --git a/matlab/ch10_ex19_ou_filtering_smoothing.m b/matlab/ch10_ex19_ou_filtering_smoothing.m new file mode 100644 index 0000000..f16b24a --- /dev/null +++ b/matlab/ch10_ex19_ou_filtering_smoothing.m @@ -0,0 +1,269 @@ +%% Examples 10.19, 10.21, 10.29, 10.33: Kalman filtering/smoothing of OU +% +% Copyright: +% 2018 - Simo Särkkä and Arno Solin +% +% License: +% This software is provided under the MIT License. See the accompanying +% LICENSE file for details. + +%% +% This is a double check of the discretization +% + dt = 0.1; + lambda = 0.3; + q = pi; + + a1 = exp(-lambda*dt); + S1 = (q/(2*lambda))*(1 - exp(-2*lambda*dt)); + + [a2,S2] = lti_disc(-lambda,1,q,dt) + +%% +% Simulate data +% + % Lock random seed + if exist('rng') % Octave doesn't have rng + rng(103,'twister'); + else + randn('state',0); + end + + dt = 0.01; + lambda = 0.5; + q = 1; + R = 0.1; + steps = 300; + + X = zeros(1,steps); + Y = zeros(1,steps); + T = zeros(1,steps); + + m0 = 0; + P0 = 1; + + x = m0 + sqrt(P0) * randn; + t = 0; + A = exp(-lambda*dt); + Q = (q - q*exp(-2*dt*lambda))/(2*lambda); + for k=1:steps + x = A*x + sqrt(Q) * randn; + if rem(k,10) == 0 + y = x + sqrt(R) * randn; + else + y = NaN; + end + t = t + dt; + X(k) = x; + Y(k) = y; + T(k) = t; + end + + figure(1); clf + plot(T,X,T,Y,'+'); + +%% +% Kalman filter (implement in 3 ways for debugging) +% + + MM = zeros(1,length(Y)); + PP = zeros(1,length(Y)); + + % This is the debug + m = m0; + P = P0; + for k=1:length(Y) + m = A*m; + P = A*P*A' + Q; + if ~isnan(Y(k)) + S = P + R; + K = P / S; + m = m + K * (Y(k) - m); + P = P - K * S * K'; + end + end + m1 = m; + P1 = P; + + % This is corresponding to Example 10.19 + m = m0; + P = P0; + eu_steps = 100; + for k=1:length(Y) + dt2 = dt / eu_steps; + for j=1:eu_steps + m = m + (-lambda * m) * dt2; + P = P + (-2 * lambda * P + q) * dt2; + end + + if ~isnan(Y(k)) + m = m + P / (P + R) * (Y(k) - m); + P = P - P^2 / (P + R); + end + MM(k) = m; + PP(k) = P; + end + m2 = m; + P2 = P; + + % This is the actual discretized one (Example 10.21) + m = m0; + P = P0; + for k=1:length(Y) + m = exp(-lambda * dt) * m; + P = exp(-2 * lambda * dt) * P + (q / (2 * lambda)) * (1 - exp(-2 * lambda * dt)); + if ~isnan(Y(k)) + m = m + P / (P + R) * (Y(k) - m); + P = P - P^2 / (P + R); + end + MM(k) = m; + PP(k) = P; + end + + [m m1 m2] + [P P1 P2] + + figure(1); + plot(T,MM,T,Y,'+',T,X,'--',T,MM-1.96*sqrt(PP),'-.',T,MM+1.96*sqrt(PP),'-.'); + +%% +% RTS smoother, again implement in 3 different ways +% + + % This is a reference solution + kf_m = MM; + kf_P = reshape(PP,[1 1 length(PP)]); + ms = kf_m(:,end); + Ps = kf_P(:,:,end); + rts_m = zeros(size(m,1),size(Y,2)); + rts_P = zeros(size(P,1),size(P,2),size(Y,2)); + rts_m(:,end) = ms; + rts_P(:,:,end) = Ps; + for k=size(kf_m,2)-1:-1:1 + mp = A*kf_m(:,k); + Pp = A*kf_P(:,:,k)*A'+Q; + Ck = kf_P(:,:,k)*A'/Pp; + ms = kf_m(:,k) + Ck*(ms - mp); + Ps = kf_P(:,:,k) + Ck*(Ps - Pp)*Ck'; + rts_m(:,k) = ms; + rts_P(:,:,k) = Ps; + end + ms1 = ms; + Ps1 = Ps; + + % This is corresponding to Example 10.29 + MMS2 = MM; + PPS2 = PP; + ms = MM(end); + Ps = PP(end); + for k=length(Y)-1:-1:1 + MMf = zeros(1,eu_steps); + PPf = zeros(1,eu_steps); + dt2 = dt / eu_steps; + m = MM(k); + P = PP(k); + for j=1:eu_steps + MMf(j) = m; + PPf(j) = P; + m = m + (-lambda * m) * dt2; + P = P + (-2 * lambda * P + q) * dt2; + end + + for j=eu_steps:-1:1 + dms = -lambda * ms + (q / PPf(j)) * (ms - MMf(j)); + dPs = 2 * (-lambda + q / PPf(j)) * Ps - q; + ms = ms - dms * dt2; + Ps = Ps - dPs * dt2; + end + + MMS2(k) = ms; + PPS2(k) = Ps; + end + ms2 = ms; + Ps2 = Ps; + + % + % This is the actual one from Example 10.33 + % + MMS = MM; + PPS = PP; + + ms = MM(end); + Ps = PP(end); + for k=length(Y)-1:-1:1 + mp = exp(-lambda * dt) * MM(k); + Pp = exp(-2 * lambda * dt) * PP(k) + (q / (2 * lambda)) * (1 - exp(-2 * lambda * dt)); + ms = MM(k) + (PP(k) * exp(-lambda * dt) / Pp) * (ms - mp); + Ps = PP(k) + (PP(k) * exp(-lambda * dt) / Pp)^2 * (Ps - Pp); + + MMS(k) = ms; + PPS(k) = Ps; + end + + [ms ms1 ms2] + + figure(1); clf; + subplot(2,1,1); + plot(T,MMS,'r-',T,rts_m,'b--',T,MMS2,'g:') + + subplot(2,1,2); + plot(T,PPS,'r-',T,squeeze(rts_P),'b--',T,PPS2,'g:') + +%% +% Plot the final figures +% + figure(1); clf; hold on + + % Plot KF mean + h1 = plot(T,MM,'-k'); + set(h1,'LineWidth',1); + set(h1,'Color',[0.5 0.5 0.5]); + + % Plot KF uncertainty + QQ1 = MM - 1.96 * sqrt(PP); + QQ2 = MM + 1.96 * sqrt(PP); + h2 = fill([T'; flipud(T')], [QQ1'; flipud(QQ2')],1); + set(h2,'EdgeColor',[.7 .7 .7],'FaceColor',[.7 .7 .7]) + + % Plot signal + h3 = plot(T,X,'-k'); + set(h3,'LineWidth',1); + + % Plot KF mean + h4 = plot(T,MM,'-k'); + set(h4,'LineWidth',1); + set(h4,'Color',[0.5 0.5 0.5]); + h5 = plot(T,Y,'+k'); + box on + xlabel('Time, $t$'), ylabel('Signal, $x(t)$') + + legend([h1, h2, h3, h5],'KF mean','KF 95\% quantiles','Signal','Observations') + + figure(2); clf; hold on + + % Plot RTS mean + h1 = plot(T,MMS,'-k'); + set(h1,'LineWidth',1); + set(h1,'Color',[0.5 0.5 0.5]); + + % Plot RTS uncertainty + QQ1 = MMS - 1.96 * sqrt(PPS); + QQ2 = MMS + 1.96 * sqrt(PPS); + h2 = fill([T'; flipud(T')], [QQ1'; flipud(QQ2')],1); + set(h2,'EdgeColor',[.7 .7 .7],'FaceColor',[.7 .7 .7]) + + % Plot signal + h3 = plot(T,X,'-k'); + set(h3,'LineWidth',1); + + % Plot RTS mean + h4 = plot(T,MMS,'-k'); + set(h4,'LineWidth',1); + set(h4,'Color',[0.5 0.5 0.5]); + + h5 = plot(T,Y,'+k'); + + box on + xlabel('Time, $t$'), ylabel('Signal, $x(t)$') + + legend([h1, h2, h3, h5],'RTS mean','RTS 95\% quantiles','Signal','Observations') diff --git a/matlab/ch11_ex05_ou_parameter_estimation.m b/matlab/ch11_ex05_ou_parameter_estimation.m new file mode 100644 index 0000000..137aade --- /dev/null +++ b/matlab/ch11_ex05_ou_parameter_estimation.m @@ -0,0 +1,381 @@ +%% Examples 11.5 and 11.9: Parameter estimation in Ornstein-Uhlenbeck +% +% Copyright: +% 2018 - Simo Särkkä and Arno Solin +% +% License: +% This software is provided under the MIT License. See the accompanying +% LICENSE file for details. + +%% +% Simulate data +% + % Lock random seed + if exist('rng') % Octave doesn't have rng + rng(1); + else + randn('state',1); + end + + dt = 0.1; + lambda = 0.5; + q = 1; + steps = 100; + + X = zeros(1,steps); + T = zeros(1,steps); + + x = 0; + t = 0; + A = exp(-lambda*dt); + Q = (q - q*exp(-2*dt*lambda))/(2*lambda); + for k=1:steps + X(k) = x; + T(k) = t; + x = A*x + sqrt(Q) * randn; + t = t + dt; + end + + figure(1); clf + plot(T,X); + +%% +% Closed-form maximum likelihood +% + lambda_est = -1/dt * log(sum(X(1:end-1) .* X(2:end)) / ... + sum(X(1:end-1) .* X(1:end-1))) + q_est = 1/(length(X)-1) * (2*lambda_est / (1 - exp(-2*lambda_est*dt))) * ... + sum((X(2:end) - exp(-lambda_est*dt) * X(1:end-1)).^2) + +%% +% Optimized maximum likelihood +% + a_f = @(lambda,q) exp(-lambda * dt); + s_f = @(lambda,q) q/(2*lambda) * (1 - exp(-2*lambda*dt)); + + ell = @(p) 0.5 * sum(log(2*pi*s_f(p(1),p(2))) + 1/s_f(p(1),p(2)) * (X(2:end) - a_f(p(1),p(2)) * X(1:end-1)).^2); + + p = fminsearch(ell,[1 1]); + lambda_est_2 = p(1) + q_est2 = p(2) + +%% +% Optimized maximum likelihood from Euler-Maruyama +% + a_f = @(lambda,q) (1 - lambda * dt); + s_f = @(lambda,q) q*dt; + + ell_em = @(p) 0.5 * sum(log(2*pi*s_f(p(1),p(2))) + 1/s_f(p(1),p(2)) * (X(2:end) - a_f(p(1),p(2)) * X(1:end-1)).^2); + + p = fminsearch(ell_em,[1 1]); + lambda_est_em = p(1) + q_est_em = p(2) + +%% +% Do some Metropolis-Hastings on the exact model. +% + % Note that we do a change of variables such that + % theta = (log q, log lambda), which + % which leads to a correction term 1/q/lambda + + % Lock random seed + if exist('rng') % Octave doesn't have rng + rng(1); + else + randn('state',2); + end + + nmcmc = 10000; + samp = []; + + theta = [0 0]; + + en = 0; + accepted = 0; + + for j=1:nmcmc + + % Draw candidate + new_theta = theta + 0.5*randn(1,2); + + % Evaluate energy + new_en = ell(exp(new_theta)) - sum(new_theta); + + % Accept or reject + if j == 1 + en = new_en; + end + a = min(1,exp(en - new_en)); + u = rand; + + if u <= a + en = new_en; + theta = new_theta; + accepted = accepted + 1; + %fprintf('Accepted: %f %f\n',theta); + else + %fprintf('Rejected\n'); + end + samp = [samp; theta]; + + if rem(j,1000) == 0 + fprintf('Acceptance rate %f\n',accepted/j); + end + end + + figure(1); clf; + subplot(1,2,1); + + plot(exp(samp(:,1)),exp(samp(:,2)),'.'); + xlabel('Parameter $\lambda$'); + ylabel('Parameter $q_c$') + + subplot(1,2,2); + [N,XE,YE] = ndhist(exp(samp(:,1:2)),20); + pcolor(XE,YE,N); + colormap(flipud(gray)) + shading flat; + xlabel('Parameter $\lambda$'); + ylabel('Parameter $q_c$') + + +%% +% Do some Metropolis-Hastings on the E-M approximation. +% + % Lock random seed + if exist('rng') % Octave doesn't have rng + rng(1); + else + randn('state',2); + end + + nmcmc = 10000; + samp_em = []; + + theta = [0 0]; + + en = 0; + accepted = 0; + + for j=1:nmcmc + + % Draw candidate + new_theta = theta + 0.5*randn(1,2); + + % Evaluate energy + new_en = ell_em(exp(new_theta)) - sum(new_theta); + + % Accept or reject + if j == 1 + en = new_en; + end + a = min(1,exp(en - new_en)); + u = rand; + + if u <= a + en = new_en; + theta = new_theta; + accepted = accepted + 1; + %fprintf('Accepted: %f %f\n',theta); + else + %fprintf('Rejected\n'); + end + samp_em = [samp_em; theta]; + + if rem(j,1000) == 0 + fprintf('Acceptance rate %f\n',accepted/j); + end + end + + figure(1); clf; + subplot(1,2,1); + plot(exp(samp_em(:,1)),exp(samp_em(:,2)),'.'); + xlabel('Parameter $\lambda$'); + ylabel('Parameter $q_c$') + + subplot(1,2,2); + [N,XE,YE] = ndhist(exp(samp_em(:,1:2)),20); + pcolor(XE,YE,N); + colormap(flipud(gray)) + shading flat; + xlabel('Parameter $\lambda$'); + ylabel('Parameter $q_c$') + + +%% +% Fokker-Planck approximation +% + L = 15; + h = 0.05; + + x_grid = (-L:h:L)'; + n = length(x_grid); + + f = -lambda * x_grid; + + Fa = q/2/h^2+f/h; + Fb = -q/h^2-f/h; + Fc = repmat(q/2/h^2,n,1); + + F = spdiags([Fa Fb Fc],-1:1,n,n); + + fprintf('Computing the matrix exponential...'); + + T_pde = expm(full(F*dt)); + pcolor(T_pde) + shading flat + + fprintf('done.\n'); + + t=0; + pp = zeros(size(x_grid)); + ind = find(x_grid == 0); + ind = find(x_grid == 0); + pp(ind) = 1/h; + for k=1:steps + pp = T_pde * pp; + pp_norm = pp ./ sum(pp) / h; + t = t + dt; + %plot(x_grid,pp_norm) + %pause(0.1); + %drawnow; + end + + figure(1); clf; + plot(x_grid,pp_norm) + + %% + % Plot the whole likelihood surfaces and MCMC results + % + [ll,qq] = meshgrid(0.01:0.01:2,0.01:0.01:2); + lh = zeros(size(ll)); + lh_em = zeros(size(ll)); + for i=1:size(ll,1) + for j=1:size(ll,2) + lh(i,j) = ell([ll(i,j) qq(i,j)]); + lh_em(i,j) = ell_em([ll(i,j) qq(i,j)]); + end + end + + figure(1); clf + subplot(2,2,1); + imagesc(ll(1,:),qq(:,1),exp(-lh)); + hold on; + + plot(lambda_est,q_est,'*') + xlabel('$\lambda$'); + ylabel('$q$'); + title('Exact likelihood'); + ax = axis; axis xy, box on + colormap(flipud(gray)) + set(gca,'Layer','Top') + set(gca,'XTick',0:.5:2.5,'YTick',0:.5:2.5) + + subplot(2,2,2); + n = size(samp,1); + for i=1:1000:n + ind = i:min(i+999,n); + plot(exp(samp(ind,1)),exp(samp(ind,2)),'.','MarkerSize',3,'Color',[.5 .5 .5]); + hold on; + end + xlabel('$\lambda$'); + ylabel('$q$'); + title('Exact MCMC'); + grid on; box on + axis(ax); + set(gca,'XTick',0:.5:2.5,'YTick',0:.5:2.5) + + subplot(2,2,3); + imagesc(ll(1,:),qq(:,1),exp(-lh_em)); + hold on; + + plot(lambda_est_em,q_est_em,'*') + xlabel('$\lambda$'); + ylabel('$q$'); + title('Pseudo-likelihood'); + ax = axis; axis xy, box on + colormap(flipud(gray)) + set(gca,'Layer','Top') + set(gca,'XTick',0:.5:2.5,'YTick',0:.5:2.5) + + subplot(2,2,4); + n = size(samp_em,1); + for i=1:1000:n + ind = i:min(i+999,n); + plot(exp(samp_em(ind,1)),exp(samp_em(ind,2)),'.','MarkerSize',3,'Color',[.5 .5 .5]); + hold on; + end + xlabel('$\lambda$'); + ylabel('$q$'); + title('Pseudo-MCMC'); + grid on; box on + axis(ax); + set(gca,'XTick',0:.5:2.5,'YTick',0:.5:2.5) + +%% +% Plot final figures +% + + % Exact likelihood + figure(1); clf; hold on; + + imagesc(ll(1,:),qq(:,1),exp(-lh)); + + h = plot(lambda_est,q_est,'*'); + set(h,'Color',[1 1 1]); + xlabel('$\lambda$'); + ylabel('$q$'); + ax = axis; axis xy, box on + colormap(flipud(gray)) + set(gca,'Layer','Top') + set(gca,'XTick',0:.5:2.5,'YTick',0:.5:2.5) + + % Exact MCMC + figure(2); clf; + n = size(samp,1); + for i=1:1000:n + ind = i:min(i+999,n); + plot(exp(samp(ind,1)),exp(samp(ind,2)),'.','MarkerSize',3,'Color',[.5 .5 .5]); + hold on; + end + xlabel('$\lambda$'); + ylabel('$q$'); + grid on; box on + axis(ax); + set(gca,'XTick',0:.5:2.5,'YTick',0:.5:2.5) + + % Pseudo-likelihood + figure(3); clf; + imagesc(ll(1,:),qq(:,1),exp(-lh_em)); + hold on; + h = plot(lambda_est_em,q_est_em,'*'); + set(h,'Color',[1 1 1]); + xlabel('$\lambda$'); + ylabel('$q$'); + ax = axis; axis xy, box on + colormap(flipud(gray)) + set(gca,'Layer','Top') + set(gca,'XTick',0:.5:2.5,'YTick',0:.5:2.5) + + % Pseudo-MCMC + figure(4); clf; + n = size(samp_em,1); + for i=1:1000:n + ind = i:min(i+999,n); + plot(exp(samp_em(ind,1)),exp(samp_em(ind,2)),'.','MarkerSize',3,'Color',[.5 .5 .5]); + hold on; + end + xlabel('$\lambda$'); + ylabel('$q$'); + grid on; box on + axis(ax); + set(gca,'XTick',0:.5:2.5,'YTick',0:.5:2.5) + + % Trajectory + figure(5); clf; + h = plot(T,X); + set(h,'Color',[0 0 0]); + set(h,'LineWidth',1); + xlabel('Time, $t$'); + ylabel('$x(t)$'); + \ No newline at end of file diff --git a/matlab/ch12_ex06_gp_regression.m b/matlab/ch12_ex06_gp_regression.m new file mode 100644 index 0000000..755d941 --- /dev/null +++ b/matlab/ch12_ex06_gp_regression.m @@ -0,0 +1,226 @@ +%% Examples 12.6 and 12.11: Batch and sequential solution to GP regression +% +% Copyright: +% 2018 - Simo Särkkä and Arno Solin +% +% License: +% This software is provided under the MIT License. See the accompanying +% LICENSE file for details. + +%% Batch GP regression + + % Lock random seed + if exist('rng') % Octave doesn't have rng + rng(0,'twister') + else + randn('state',0); + end + + % Define parameters + sigma2 = 0.1^2; + magnSigma2 = 1^2; + ell = .1; + + % Define covariance function: Matern 3/2 + C = @(x,y) magnSigma2*(1+sqrt(3)*abs(x(:)-y(:)')/ell).* ... + exp(-sqrt(3)*abs(x(:)-y(:)')/ell); + + % Number of inputs + n = 10; + to = sort(rand(n,1)); + + % Simulate data + Koo = C(to,to); + cK = chol(Koo); + yo = cK'*randn(size(cK,1),1) + sqrt(sigma2)*randn(n,1); + + % Evaluate test points + t = linspace(0,1,500)'; + Kto = C(t,to); + Ktt = C(t,t); + + % Solve the GP regression problem + L = chol(Koo+sigma2*eye(n),'lower'); + m_gp = Kto*((Koo+sigma2*eye(n))\yo); + cov_gp = Ktt - Kto*((Koo+sigma2*eye(n))\Kto'); + var_gp = diag(cov_gp); + + % Samples form the posterior GP + ns = 5; + cK = chol(cov_gp); + xs = m_gp + cK'*randn(size(cK,1),ns); + + % Plot + figure(1); clf; hold on + + % Plot mean + h1 = plot(t,m_gp,'-k'); + + % Plot uncetainty + h2 = fill([t; flipud(t)],[m_gp+1.96*sqrt(var_gp); flipud(m_gp-1.96*sqrt(var_gp))],1); + set(h2,'EdgeColor',[.7 .7 .7],'FaceColor',[.7 .7 .7]) + + % Plot sample trajectories + h4 = plot(t,xs,'-','Color',[.5 .5 .5],'LineWidth',.5); + + % Plot mean + plot(t,m_gp,'-k'); + + % Plot observations + h3 = plot(to,yo,'+k'); + + % Axis options + set(gca,'layer','top') + box on + ylim([-3 4]) + xlabel('Input, $\xi$'), ylabel('Output, $x(\xi)$') + legend([h1, h2, h3, h4(1)],'Mean','95\% quantiles','Observations','Samples') + + +%% GP regression by Kalman filtering and RTS smoothing + + % Define the LTI SDE model: Matern 3/2 + lambda = sqrt(3)/ell; + F = [0 1; -lambda^2 -2*lambda]; + L = [0; 1]; + Q = 4*lambda^3*magnSigma2; + H = [1 0]; + Pinf = magnSigma2*diag([1, lambda^2]); + + % Combine test and observation times + [tt,ind] = sort([to; t]); + + % Time steps with observations + doUpdate = (ind<=n); + + % Set initial mean and covariance + m = zeros(size(F,1),1); + P = Pinf; + + % Allocate space for results + MF = zeros(size(m,1),numel(tt)); + PF = zeros(size(m,1),size(m,1),numel(tt)); + MP = zeros(size(m,1),numel(tt)); + PP = zeros(size(m,1),size(m,1),numel(tt)); + GS = zeros(size(m,1),size(m,1),numel(tt)); + + % Run the Kalman filter + for k=1:numel(tt) + + % Time step + if k==1 + dt = 0; + else + dt = tt(k)-tt(k-1); + end + + % Solve the LTI SDE for this step + A = expm(F*dt); + Q = Pinf - A*Pinf*A'; + + % Kalman prediction + mp = A*m; + Pp = A*P*A' + Q; + + % Pre-calculate smoother gain + Gs = (P*A')/Pp; + + % Do Kalman update if there is an observation on this step + if doUpdate(k) + v = yo(ind(k)) - H*m; + S = H*P*H' + sigma2; + K = P*H'/S; + m = m + K*v; + P = P - K*H*P; + else + m = mp; + P = Pp; + end + + % Store results + MF(:,k) = m; + PF(:,:,k) = P; + MP(:,k) = mp; + PP(:,:,k) = Pp; + GS(:,:,k) = Gs; + + end + + % Run the RTS smoother + MS = MF; + PS = PF; + ms = MS(:,end); + Ps = PS(:,:,end); + for j=numel(tt)-1:-1:1 + ms = MF(:,j) + GS(:,:,j+1)*(ms-MP(:,j+1)); + Ps = PF(:,:,j) + GS(:,:,j+1)*(Ps - PP(:,:,j+1))*GS(:,:,j+1)'; + MS(:,j) = ms; + PS(:,:,j) = Ps; + end + + % Extract values + m_f = (H*MF)'; + var_f = arrayfun(@(k) H*PF(:,:,k)*H',1:numel(tt))'; + m_s = (H*MS)'; + var_s = arrayfun(@(k) H*PS(:,:,k)*H',1:numel(tt))'; + + % Only the test inputs + ii = ind>n; + m_f = m_f(ii); + var_f = var_f(ii); + m_s = m_s(ii); + var_s = var_s(ii); + + + % Plot + figure(2); clf; hold on + + % Plot uncetainty + h2 = fill([t; flipud(t)], ... + [m_f+1.96*sqrt(var_f); flipud(m_f-1.96*sqrt(var_f))],1); + set(h2,'EdgeColor',[.7 .7 .7],'FaceColor',[.7 .7 .7]) + + plot(t,m_gp+1.96*sqrt(var_gp),'--k', ... + t,m_gp-1.96*sqrt(var_gp),'--k'); + + % Plot mean + plot(t,m_f,'-k') + + % Plot observations + plot(to,yo,'+k') + + % Axis options + set(gca,'layer','top') + box on + xlabel('Input, $t$'), ylabel('Output, $x(t)$') + ylim([-3 4]) + + % Plot + figure(3); clf; hold on + + % Plot mean + h1 = plot(t,m_s,'-k'); + + % Plot uncetainty + h2 = fill([t; flipud(t)], ... + [m_s+1.96*sqrt(abs(var_s)); flipud(m_s-1.96*sqrt(abs(var_s)))],1); + set(h2,'EdgeColor',[.7 .7 .7],'FaceColor',[.7 .7 .7]) + + % Plot mean + plot(t,m_s,'-k') + + % Plot observations + h3 = plot(to,yo,'+k'); + + h4 = plot(t,m_gp+1.96*sqrt(var_gp),'--k', ... + t,m_gp-1.96*sqrt(var_gp),'--k'); + + % Axis options + set(gca,'layer','top') + box on + xlabel('Input, $t$'), ylabel('Output, $x(t)$') + ylim([-3 4]) + legend([h1, h2, h3, h4(1)],'Mean','95\% quantiles','Observations','Batch solution') + + + \ No newline at end of file diff --git a/matlab/ch12_ex12_gp_approximation_of_drift.m b/matlab/ch12_ex12_gp_approximation_of_drift.m new file mode 100644 index 0000000..33592b8 --- /dev/null +++ b/matlab/ch12_ex12_gp_approximation_of_drift.m @@ -0,0 +1,111 @@ +%% Example 12.12: GP approximation of drift functions (double-well) +% +% Copyright: +% 2018 - Simo Särkkä and Arno Solin +% +% License: +% This software is provided under the MIT License. See the accompanying +% LICENSE file for details. + +%% + + % Lock random seed + if exist('rng') % Octave doesn't have rng + rng(2,'twister') + else + randn('state',2); + end + + % The double-well model + f = @(x,t) 4*(x-x.^3); + L = @(x,t) 1; + q = 1; + x0 = 0; + + dt = 0.01; + t = 0:dt:(20000-1)*dt; + + % Simulate + [x,tspan] = eulermaruyama(f,L,t,x0,q); + + % Observed + ty = t(1:20:end); + y = x(1:20:end); + + %ty = ty(1:1000); + %y = y(1:1000); + Dt = ty(2)-ty(1); + + % Show + figure(1); clf; hold on + plot(t,x,'-','Color',[.5 .5 .5]) + plot(ty,y,'-k') + + % Plot the true f + figure(2); clf + subplot(211) + z = linspace(-2.5,2.5,100); + plot(z,f(z,0)) + xlim([-2.5 2.5]); ylim([-5 5]) + subplot(212) + hist(y,z) + xlim([-2.5 2.5]); + + % Define GP model + sigma_lin2 = 10^2; + sigma_se2 = 20^2; + ell = 1; + k = @(x,y) sigma_lin2*x(:)*y(:)' + sigma_se2*exp(-(x(:)-y(:)').^2/2/ell^2); + + % GP prediction + d = (y(2:end)-y(1:end-1))'/Dt; + gpf_mean = @(x) k(x,y(1:end-1)) * ((k(y(1:end-1),y(1:end-1))+q/Dt*eye(numel(y)-1))\d); + gpf_cov = @(x) k(x,x) - k(x,y(1:end-1))*((k(y(1:end-1),y(1:end-1))+q/Dt*eye(numel(y)-1))\k(y(1:end-1),x)); + + figure(3); clf; hold on + z = linspace(-2.5,2.5,100); + + % Predict + m = gpf_mean(z); + v = diag(gpf_cov(z)); + + % Plot mean + h1 = plot(z,m,'-k'); + + % Plot uncertainty + h2 = fill([z'; flipud(z')],[m+1.96*sqrt(v); flipud(m-1.96*sqrt(v))],1); + set(h2,'EdgeColor',[.7 .7 .7],'FaceColor',[.7 .7 .7]) + + % Plot mean + plot(z,m,'-k'); + + % Plot ground-truth + h3 = plot(z,f(z,0),'--k'); + + % Limits + box on + xlim([-2.5 2.5]); ylim([-5 5]) + xlabel('Input, $x$') + ylabel('$f(x)$') + + % Ticks + set(gca,'Layer','top') + set(gca,'XTick',-2:2,'YTick',-4:2:4) + + % Legend + legend([h1 h2 h3],'Mean','95\% quantiles','Ground-truth drift') + + figure(4); clf + plot(ty,y,'-k') + + % Limits + box on + ylim([-2.5 2.5]); + xlabel('Time, $t$') + ylabel('$x(t)$') + + % Ticks + set(gca,'Layer','top') + set(gca,'YTick',-2:2) + + \ No newline at end of file diff --git a/matlab/euler.m b/matlab/euler.m new file mode 100644 index 0000000..f284c09 --- /dev/null +++ b/matlab/euler.m @@ -0,0 +1,49 @@ +function [x,tspan] = euler(f,tspan,x0) +%% EULER - Numerical ODE solver: Euler's method +% +% Syntax: +% [x,tspan] = euler(f,tspan,x0) +% +% In: +% f - Function handle, f(x,t) +% tspan - Time steps to simulate, [t0,...,tend] +% x0 - Initial condition +% +% Out: +% x - Solved values +% tspan - Time steps +% +% Description: +% Integrates the system of differential equations +% x' = f(x,t), for x(0) = x0 +% over the time interval defined in tspan. +% +% Copyright: +% 2018 - Simo Särkkä and Arno Solin +% +% License: +% This software is provided under the MIT License. See the accompanying +% LICENSE file for details. + +%% + + % Number of steps + steps = numel(tspan); + + % Allocate space + x = zeros(size(x0,1),steps); + + % Initial state + x(:,1) = x0; + + % Iterate + for k=2:steps + + % Time discretization + dt = tspan(k)-tspan(k-1); + + % Step + x(:,k) = x(:,k-1) + f(x(:,k-1),tspan(k-1))*dt; + + end + \ No newline at end of file diff --git a/matlab/eulermaruyama.m b/matlab/eulermaruyama.m new file mode 100644 index 0000000..e068981 --- /dev/null +++ b/matlab/eulermaruyama.m @@ -0,0 +1,63 @@ +function [x,tspan] = eulermaruyama(f,L,tspan,x0,Q) +%% EULERMARUYAMA - Numerical SDE solver: The Euler-Maruyama method +% +% Syntax: +% [x,tspan] = eulermaruyama(f,L,tspan,x0,Q) +% +% In: +% f - Drift function, f(x,t) +% L - Diffusion function, L(x,t) +% tspan - Time steps to simulate, [t0,...,tend] +% x0 - Initial condition +% Q - Spectral density (default: standard Brownian motion) +% +% Out: +% x - Solved values +% tspan - Time steps +% +% Description: +% Integrates the system of stochatic differential equations +% dx = f(x,t) dt + L(x,t) dbeta, for x(0) = x0 +% over the time interval defined in tspan. +% +% Copyright: +% 2018 - Simo Särkkä and Arno Solin +% +% License: +% This software is provided under the MIT License. See the accompanying +% LICENSE file for details. + +%% + + % Check if Q given + if nargin<5 || isempty(Q), Q = eye(size(L(x0,tspan(1)),2)); end + + % Cholesky factor of Q + cQ = chol(Q,'lower'); + + % Number of steps + steps = numel(tspan); + + % Allocate space + x = zeros(size(x0,1),steps); + + % Initial state + x(:,1) = x0; + + % Iterate + for k=2:steps + + % Time discretization + dt = tspan(k)-tspan(k-1); + + % Increment + db = sqrt(dt)*cQ*randn(size(Q,1),1); + + % Step + x(:,k) = x(:,k-1) + ... + f(x(:,k-1),tspan(k-1))*dt + ... + L(x(:,k-1),tspan(k-1))*db; + + end + + diff --git a/matlab/eulermaruyama_weak.m b/matlab/eulermaruyama_weak.m new file mode 100644 index 0000000..5b440a5 --- /dev/null +++ b/matlab/eulermaruyama_weak.m @@ -0,0 +1,74 @@ +function [x,tspan] = eulermaruyama_weak(f,L,tspan,x0,Q,gaussian) +%% EULERMARUYAMA_WEAK - Numerical SDE solver: The weak Euler-Maruyama method +% +% Syntax: +% [x,tspan] = eulermaruyama_weak(f,L,tspan,x0,Q,gaussian) +% +% In: +% f - Drift function, f(x,t) +% L - Diffusion function, L(x,t) +% tspan - Time steps to simulate, [t0,...,tend] +% x0 - Initial condition +% Q - Spectral density (default: standard Brownian motion) +% gaussian - Use Gaussian increments (default: false) +% +% Out: +% x - Solved values +% tspan - Time steps +% +% Description: +% Integrates the system of stochatic differential equations +% dx = f(x,t) dt + L(x,t) dbeta, for x(0) = x0 +% over the time interval defined in tspan. +% +% Copyright: +% 2018 - Simo Särkkä and Arno Solin +% +% License: +% This software is provided under the MIT License. See the accompanying +% LICENSE file for details. + +%% + + % Check if Q given + if nargin<5 || isempty(Q), Q = eye(size(L(x0,tspan(1)),2)); end + + % Check the sort of increments to use + if nargin<6 || isempty(gaussian), gaussian = false; end + + % Cholesky factor of Q + cQ = chol(Q,'lower'); + + % Number of steps + steps = numel(tspan); + + % Allocate space + x = zeros(size(x0,1),steps); + + % Initial state + x(:,1) = x0; + + % Pre-calculate random numbers + if gaussian + R = randn(size(Q,1),steps); + else + R = (rand(size(Q,1),steps)>.5)*2-1; + end + + % Iterate + for k=2:steps + + % Time discretization + dt = tspan(k)-tspan(k-1); + + % Increment + db = sqrt(dt)*cQ*R(:,k); + + % Step + x(:,k) = x(:,k-1) + ... + f(x(:,k-1),tspan(k-1))*dt + ... + L(x(:,k-1),tspan(k-1))*db; + + end + + diff --git a/matlab/herm_exp.m b/matlab/herm_exp.m new file mode 100644 index 0000000..604ad9d --- /dev/null +++ b/matlab/herm_exp.m @@ -0,0 +1,53 @@ +function [ herm_x, ser, c_list, H_list ] = herm_exp(f, df, d2f, d3f, d4f, d5f, dt, x, x0) +%% HERM_EXP - Hermite expansion +% +% Syntax: +% [herm_x, ser, c_list, H_list] = herm_exp(f, df, d2f, d3f, d4f, d5f, dt, x, x0) +% +% Description: +% Supporting file for 'ch09_ex18_hermite_expansion.m'. +% +% Copyright: +% 2018 - Simo Särkkä and Arno Solin +% +% License: +% This software is provided under the MIT License. See the accompanying +% LICENSE file for details. + +%% + + c_list = {}; + c_list{end+1} = ones(size(x)); + c_list{end+1} = -f.*dt^(1/2) - (2.*f.*df + d2f) .* dt^(3/2) / 4 ... + - (4.*f.*df.^2 + 4.*f.^2.*d2f + 6.*df.*d2f + 4.*f.*d3f + d4f) .* dt^(5/2) / 24; + c_list{end+1} = (f.^2 + df).*dt/2 + (6.*f.^2.*df + 4.*df.^2 + 7.*f.*d2f + 2.*d3f).*dt^2/12 ... + + (28 .* f.^2 .* df.^2 + 28 .* f.^2 .* d3f + 16 .* df.^3 + 16 .* f.^3 .* d2f + 88.*f.*df.*d2f ... + + 21.*d2f.^2 + 32.*df.*d3f + 16.*f.*d4f + 3.*d5f) .* dt^3 / 96; + c_list{end+1} = -(f.^3 + 3.*f.*df + d2f).*dt^(3/2)/6 - (12.*f.^3.*df + 28.*f.*df.^2 ... + + 22.*f.^2.*d2f + 24.*df.*d2f + 14.*f.*d3f + 3.*d4f) .* dt^(5/2) / 48; + c_list{end+1} = (f.^4 + 6.*f.^2.*df + 3.*df.^2 + 4.*f.*d2f + d3f).*dt^2/24 ... + + (20.*f.^4.*df + 50.*f.^3.*d2f + 100.*f.^2.*df.^2 + 50.*f.^2.*d3f + 23.*f.*d4f ... + + 180.*f.*df.*d2f + 40.*df.^3 + 34.*d2f.^2 + 52.*df.*d3f + 4.*d5f).*dt^3/240; + c_list{end+1} = -(f.^5 + 10.*f.^3.*df + 15.*f.*df.^2 + 10.*f.^2.*d2f ... + + 10.*df.*d2f + 5.*f.*d3f + d4f).*dt^(5/2)/120; + c_list{end+1} = (f.^6 + 15.*f.^4.*df + 15.*df.^3 + 20.*f.^3.*d2f + 15.*df.*d3f + 45.*f.^2.*df.^2 ... + + 10.*d2f.^2 + 15.*f.^2.*d3f + 60.*f.*df.*d2f + 6.*f.*d4f + d5f).*dt^3/720; + + H_list = {}; + H_list{end+1} = ones(size(x)); + H_list{end+1} = -((x-x0)/sqrt(dt)); + H_list{end+1} = -1 + ((x-x0)/sqrt(dt)).^2; + H_list{end+1} = 3.*((x-x0)/sqrt(dt)) - ((x-x0)/sqrt(dt)).^3; + H_list{end+1} = 3 - 6.*((x-x0)/sqrt(dt)).^2 + ((x-x0)/sqrt(dt)).^4; + H_list{end+1} = -15.*((x-x0)/sqrt(dt)) + 10.*((x-x0)/sqrt(dt)).^3 - ((x-x0)/sqrt(dt)).^5; + H_list{end+1} = -15 + 45.*((x-x0)/sqrt(dt)).^2 - 15.*((x-x0)/sqrt(dt)).^4 + ((x-x0)/sqrt(dt)).^6; + + ser = zeros(size(x)); + for i=1:length(c_list) + ser = ser + c_list{i} .* H_list{i}; + end + + herm_x = exp(-(x-x0).^2/(2*dt)) / sqrt(2*pi*dt) .* ser; + +end + diff --git a/matlab/heun.m b/matlab/heun.m new file mode 100644 index 0000000..e2cf0f7 --- /dev/null +++ b/matlab/heun.m @@ -0,0 +1,52 @@ +function [x,tspan] = heun(f,tspan,x0) +%% HEUN - Numerical ODE solver: Heun's method +% +% Syntax: +% [x,tspan] = heun(f,tspan,x0) +% +% In: +% f - Function handle, f(x,t) +% tspan - Time steps to simulate, [t0,...,tend] +% x0 - Initial condition +% +% Out: +% x - Solved values +% tspan - Time steps +% +% Description: +% Integrates the system of differential equations +% x' = f(x,t), for x(0) = x0 +% over the time interval defined in tspan. +% +% Copyright: +% 2018 - Simo Särkkä and Arno Solin +% +% License: +% This software is provided under the MIT License. See the accompanying +% LICENSE file for details. + +%% + + % Number of steps + steps = numel(tspan); + + % Allocate space + x = zeros(size(x0,1),steps); + + % Initial state + x(:,1) = x0; + + % Iterate + for k=2:steps + + % Time discretization + dt = tspan(k)-tspan(k-1); + + % Supporting values + fx = f(x(:,k-1),tspan(k-1)); + tx2 = x(:,k-1) + fx*dt; + + % Step + x(:,k) = x(:,k-1) + (fx + f(tx2,tspan(k-1)+dt))/2*dt; + + end \ No newline at end of file diff --git a/matlab/impliciteuler.m b/matlab/impliciteuler.m new file mode 100644 index 0000000..0db3479 --- /dev/null +++ b/matlab/impliciteuler.m @@ -0,0 +1,51 @@ +function [x,tspan] = impliciteuler(f,tspan,x0) +%% IMPLICITEULER - Numerical ODE solver: Implicit Euler method +% +% Syntax: +% [x,tspan] = euler(f,tspan,x0) +% +% In: +% f - Function handle, f(x,t) +% tspan - Time steps to simulate, [t0,...,tend] +% x0 - Initial condition +% +% Out: +% x - Solved values +% tspan - Time steps +% +% Description: +% Integrates the system of differential equations +% x' = f(x,t), for x(0) = x0 +% over the time interval defined in tspan. +% +% Copyright: +% 2018 - Simo Särkkä and Arno Solin +% +% License: +% This software is provided under the MIT License. See the accompanying +% LICENSE file for details. + +%% + + % Number of steps + steps = numel(tspan); + + % Allocate space + x = zeros(size(x0,1),steps); + + % Initial state + x(:,1) = x0; + + % Iterate + for k=2:steps + + % Time discretization + dt = tspan(k)-tspan(k-1); + + % Step solve the algebraic problem + x(:,k) = fsolve(@(z) x(:,k-1) + f(z,tspan(k))*dt - z, ... + x(:,k-1), optimset('display','none')); + + end + + \ No newline at end of file diff --git a/matlab/leapfrog.m b/matlab/leapfrog.m new file mode 100644 index 0000000..a58c12d --- /dev/null +++ b/matlab/leapfrog.m @@ -0,0 +1,73 @@ +function [x,tspan] = leapfrog(f,s,eta,tspan,x0,Q) +%% LEAPFROG - Numerical SDE solver: The simple leapfrog Verlet method +% +% Syntax: +% [x,tspan] = leapfrog(f,L,tspan,x0,Q) +% +% In: +% f - Function f(x), see below +% s - Function s(x), see below +% nu - Scalar constant, see below +% tspan - Time steps to simulate, [t0,...,tend] +% x0 - Initial condition +% Q - Spectral density (default: standard Brownian motion) +% +% Out: +% x - Solved values +% tspan - Time steps +% +% Description: +% Integrates the second-order stochatic differential equation +% \ddot{x} = f(x) - \nu s^2(x) \dot{x} + s(x) w(t), for x(0) = x0 +% over the time interval defined in tspan. +% +% References: +% [1] Burrage, Lenane, and Lythe (2007). NUMERICAL METHODS FOR +% SECOND-ORDER STOCHASTIC DIFFERENTIAL EQUATIONS. SIAM J. +% SCI. COMPUT. +% +% Copyright: +% 2018 - Simo Särkkä and Arno Solin +% +% License: +% This software is provided under the MIT License. See the accompanying +% LICENSE file for details. + +%% + + % Check if Q given + if nargin<6 || isempty(Q), Q = 1; end + + % Square-root of Q + cQ = sqrt(Q); + + % Number of steps + steps = numel(tspan); + + % Allocate space + x = zeros(size(x0,1),steps); + + % Initial state + x(:,1) = x0; + + % Iterate + for k=1:steps-1 + + % Time discretization + dt = tspan(k+1)-tspan(k); + + % Gaussian increment + db = sqrt(dt)*cQ*randn(1); + + % Supporting half-step value + hx = x(1,k) + 1/2*x(2,k)*dt; + + % Velocity + x(2,k+1) = x(2,k) - eta*s(hx)^2*x(2,k)*dt + f(hx)*dt + s(hx)*db; + + % Position + x(1,k+1) = hx + 1/2*x(2,k+1)*dt; + + end + + \ No newline at end of file diff --git a/matlab/lti_disc.m b/matlab/lti_disc.m new file mode 100644 index 0000000..8416372 --- /dev/null +++ b/matlab/lti_disc.m @@ -0,0 +1,67 @@ +function [A,Q] = lti_disc(F,L,Q,dt) +% LTI_DISC Equivalent discrete-time solution of an LTI SDE +% +% Syntax: +% [A,Q] = lti_disc(F,L,Qc,dt) +% +% In: +% F - NxN Feedback matrix +% L - NxL Noise effect matrix (optional, default identity) +% Qc - LxL Diagonal Spectral Density (optional, default zeros) +% dt - Time step (optional, default 1) +% +% Out: +% A - Transition matrix +% Q - Discrete process covariance matrix +% +% Description: +% Equivalent discrete-time solution of an LTI SDE of form +% +% dx/dt = F x + L w, +% +% where w(t) is a white-noise process with spectral density Qc. +% Results in the following model: +% +% x[k] = A x[k-1] + q, q ~ N(0,Q). +% +% Can be used for integrating the model exactly over time steps, +% which are multiples of dt. +% +% Copyright: +% 2019 - Simo Särkkä and Arno Solin +% +% License: +% This software is provided under the MIT License. See the accompanying +% LICENSE file for details. + + % Check number of arguments and defaults + if nargin < 1 + error('Too few arguments'); + end + if nargin < 2 + L = []; + end + if nargin < 3 + Q = []; + end + if nargin < 4 + dt = []; + end + if isempty(L) + L = eye(size(F,1)); + end + if isempty(Q) + Q = zeros(size(F,1),size(F,1)); + end + if isempty(dt) + dt = 1; + end + + % Closed-form integration of the transition matrix + A = expm(F*dt); + + % Closed-form integration of the covariance matrix + n = size(F,1); + Phi = [F L*Q*L'; zeros(n,n) -F']; + AB = expm(Phi*dt)*[zeros(n,n);eye(n)]; + Q = AB(1:n,:)*A'; % A' = inv(AB((n+1):(2*n),:)); \ No newline at end of file diff --git a/matlab/ndhist.m b/matlab/ndhist.m new file mode 100644 index 0000000..37fb086 --- /dev/null +++ b/matlab/ndhist.m @@ -0,0 +1,173 @@ +function [H,varargout] = ndhist(X,bins,mins,maxs,cdim) +%NDHIST Normalized histogram of N-dimensional data +% +% [H,P] = ndhist(X,[bins,mins,maxs]) +% [H,x1,x2,...] = ndhist(X,[bins,mins,maxs]) +% +% Returns normalized N-dimensional histogram H and the bin center +% coordinates in P or variables x1,x2,... which all are of the +% same size as H and contain the coordinates in same form as +% output of ndgrid. +% +% Histogram is calculated for the points in argument X. +% bins is a vector containing number of bins for each dimension. +% Optional arguments mins and maxx specify the maximum allowed +% values for each dimension. +% +% If no output arguments are given the function plots +% graph of the histogram. +% +% Examples: +% >> X = gamrnd(3,3,1000,1); +% >> [H,x] = ndhist(X,50); +% >> bar(x,H); +% +% >> X = mvnrnd([1 1],eye(2),1000) + round(rand(1000,1))*[-2 -2]; +% >> [H,x,y] = ndhist(X,[20 20]); +% >> surf(x,y,H); +% +% Copyright: +% 2019 - Simo Särkkä and Arno Solin +% +% License: +% This software is provided under the MIT License. See the accompanying +% LICENSE file for details. + +% Check input and output arguments +if nargin < 2 + bins = []; +end +if nargin < 3 + mins = []; +end +if nargin < 4 + maxs = []; +end +if nargin < 5 + cdim = 1; + if size(X,1)==1 + X = X'; + end +end + +if isempty(bins) + bins = 10; +end +if isempty(mins) + mins = min(X); +end +if isempty(maxs) + maxs = max(X); +end + +tmp = (maxs - mins) ./ bins; +mins = mins - tmp/2; +maxs = maxs + tmp/2; +dx = (maxs - mins) ./ bins; + +if (nargout~=0) && (nargout~=1) && (nargout~=2) && (nargout~=size(X,2)+1) + error('Illegal number of output arguments'); +end + +bins = reshape(bins,1,prod(size(bins))); +if size(bins,2)~=size(X,2) + if size(bins,2)==1 + bins = bins * ones(1,size(X,2)); + else + error('Wrong number of bin sizes'); + end +end + +% Create histogram for this dimension and +% recursively for each sub-dimension. +[tmpx,tmpi] = sort(X(:,cdim)); +X = X(tmpi,:); +if cdim < size(X,2) + H = zeros(bins(cdim),prod(bins(cdim+1:end))); +else + H = zeros(bins(cdim),1); +end +bnd = mins(cdim)+dx(cdim); +beg = 1; +ind = 1; +k = 1; + +while (k <= size(X,1)) && (X(k,cdim) < mins(cdim)) + k = k + 1; +end +beg = k; +while (k <= size(X,1)) && (X(k,cdim) <= maxs(cdim)) && (ind <= size(H,1)) + if X(k,cdim) <= bnd + k = k + 1; + else + if k-beg > 0 + if cdim < size(X,2) + H(ind,:) = ndhist(X(beg:k-1,:),bins,mins,maxs,cdim+1); + else + H(ind) = k-beg; + end + end + bnd = bnd + dx(cdim); + ind = ind + 1; + beg = k; + end +end + +% Handle the last bin +if (ind <= size(H,1)) + if cdim < size(X,2) + H(ind,:) = ndhist(X(beg:k-1,:),bins,mins,maxs,cdim+1); + else + H(ind) = k-beg; + end +end + +% Reshape H to linear or multidimensional and +% determine point positions and set output variables +if cdim == 1 + if size(X,2)==1 + H = H(:); + x = (0:(bins-1))'*(maxs-mins)/bins+mins+dx/2; + varargout{1} = x; + else + H = reshape(H,bins); + k = 0; + tmp = cell(size(X,2),1); + for i=1:size(X,2) + x = (0:(bins(i)-1))'/bins(i)*(maxs(i)-mins(i))+mins(i)+dx(i)/2; + args = cell(size(X,2),1); + for j=1:size(X,2) + args{j} = x; + end + tmp{i} = shiftdim(ndgrid(args{:}),k); + k = mod(k + size(X,2)-1,size(X,2)); + end + if nargout==2 + P = zeros(size(X,2),prod(bins)); + for i=1:size(X,2) + P(i,:) = reshape(tmp{i},1,prod(bins)); + end + varargout{1} = reshape(P,[size(X,2) bins]); + else + varargout = tmp; + end + end + c = sum(reshape(H,1,prod(bins))); + H = H/c; +else + H = reshape(H,1,prod(bins(cdim:end))); +end + +% If there were no output arguments, draw a graph +if nargout==0 + if size(X,2)==1 + bar(varargout{1},H,'hist'); + elseif size(X,2)==2 + surf(varargout{1},varargout{2},H); + else + error('Unable to draw >2 dimensional densities'); + end + clear H; + clear varargout; +end + diff --git a/matlab/newquiver.m b/matlab/newquiver.m new file mode 100644 index 0000000..aee1dc4 --- /dev/null +++ b/matlab/newquiver.m @@ -0,0 +1,123 @@ +function newquiver(X,Y,U,V,varargin) +% NEWQUIVER - Plot velocity vectors as arrows +% +% Syntax: +% NEWQUIVER(X,Y,U,V,options) +% +% In: +% X - Grid of x coordinates +% Y - Grid of y coordinates +% U - Vector direction (x component) +% V - Vector direction (y component) +% options - name-value pairs (see below) +% +% Out: +% (none) +% +% Description: +% NEWQUIVER(X,Y,U,V) plots velocity vectors as arrows with components +% (u,v) at the points (x,y). The matrices X,Y,U,V must all be the same +% size and contain corresponding position and velocity components +% (X and Y can also be vectors to specify a uniform grid). +% NEWQUIVER automatically scales the arrows to fit within the grid. +% +% Additional option-value pairs can be given. +% Accepted options/values are: +% * scale : A two-element vector for scaling the shape +% * color : Fill color (rgbkcy...) (default 'k': black) +% * EdgeColor : Edge color of shape (default 'k': black) +% * LineWidth : Line width (default 0.5) +% * X : Arrow shape as a two-column vector (default: arrow) +% +% See also: +% quiver +% +% Copyright: +% 2012-2018 - Simo Särkkä and Arno Solin +% +% License: +% This software is provided under the MIT License. See the accompanying +% LICENSE file for details. + +%% Define arrow shape + + % Define arrow shape + arrow.X = [8.7185878, 4.0337352 % Upper right corner + -2.2072895, 0 % Arrow tip + 8.7185878, -4.0337352 % Lower right corner + 6.9831476, -1.6157441 % Inner lower corner + 6.9831476, -0.8 % Inner lower shaft start + 28, -0.8 % Lower shaft end + 28, 0.8 % Upper shaft end + 6.9831476, 0.8 % Inner upper shaft start + 6.9831476, 1.6157441 % Inner upper corner + 8.7185878, 4.0337352]; % Upper right corner + + % Flip arrow left to right + arrow.X(:,1) = -arrow.X(:,1); + + % Normalize length to one and center + arrow.X(:,1) = arrow.X(:,1)-min(arrow.X(:,1)); + arrow.X = arrow.X/max(arrow.X(:,1)); + arrow.X(:,1) = arrow.X(:,1)-.5; + + +%% Set default options + + % Default options structure + defaultopt.scale = .5*[min(X(X(:)>min(X(:))))-min(X(:)) ... + min(Y(Y(:)>min(Y(:))))-min(Y(:))]; + defaultopt.color = 'k'; + defaultopt.EdgeColor = 'k'; + defaultopt.LineWidth = .5; + defaultopt.X = arrow.X; + + +%% Get options + + % If there are options submitted + if length(varargin) > 1 + + % Options + options = cell2struct(varargin(2:2:end),varargin(1:2:end),2); + + % Merge the default and submitted options + fnames = fieldnames(defaultopt); + for i=1:length(fnames) + %val = getfield(options,fnames{i}); + if ~isfield(options,fnames{i}) + options=setfield(options,fnames{i},getfield(defaultopt,fnames{i})); + end + end + else + options = defaultopt; + end + + % Make sure there is a scaling factor in both directions + if numel(options.scale)~=2, options.scale = [1 1]*options.scale(1); end + + +%% Get arrow directions + + % Calculate arrow direction + TH = atan2(V,U); + + hold on + + for i=1:numel(X) + + % Roatate + Z = [cos(TH(i)) -sin(TH(i)); sin(TH(i)) cos(TH(i))]*options.X'; + + % Scale + Z = bsxfun(@times,Z,options.scale(:)); + + % Move + Z = bsxfun(@plus,Z,[X(i); Y(i)]); + + % Plot + patch(Z(1,:),Z(2,:),1,'FaceColor',options.color, ... + 'EdgeColor',options.EdgeColor, ... + 'LineWidth',options.LineWidth,'LineStyle','-') + + end diff --git a/matlab/rejectionsample.m b/matlab/rejectionsample.m new file mode 100644 index 0000000..b648580 --- /dev/null +++ b/matlab/rejectionsample.m @@ -0,0 +1,54 @@ +function x = rejectionsample(q1,epsilon) +%% REJECTIONSAMPLE - Rejection sampling with Gaussian proposals +% +% Syntax: +% x = rejectionsample(q1,epsilon) +% +% In: +% q1 - Function handle +% epsilon - Bounding parameter +% +% Out: +% x - Sample +% +% Description: +% Rejection sampling with Gaussian proposals. +% +% Copyright: +% 2018 - Simo Särkkä and Arno Solin +% +% License: +% This software is provided under the MIT License. See the accompanying +% LICENSE file for details. + +%% Set defaults + + % Scale constant + if nargin<2 || isempty(epsilon) + epsilon = 1; + end + + % Proposal distribution (standard Normal) + % q2 = @(x) normpdf(x,0,1); + q2 = @(x) exp(-0.5 * x.^2) / sqrt(2*pi); + + +%% Draw a sample + + while true + + % Sample x from q2 + x = randn(1); + + % Sample a uniform random variable + u = rand(1); + + % Accept or reject x + if u < epsilon*q1(x)/q2(x) + return + end + + end + + + diff --git a/matlab/rk4simple.m b/matlab/rk4simple.m new file mode 100644 index 0000000..d87382c --- /dev/null +++ b/matlab/rk4simple.m @@ -0,0 +1,56 @@ +function [x,tspan] = rk4simple(f,tspan,x0) +%% RK4SIMPLE - Numerical ODE solver: The fourth-order Runge-Kutta method +% +% Syntax: +% [x,tspan] = rk4simple(f,tspan,x0) +% +% In: +% f - Function handle, f(x,t) +% tspan - Time steps to simulate, [t0,...,tend] +% x0 - Initial condition +% +% Out: +% x - Solved values +% tspan - Time steps +% +% Description: +% Integrates the system of differential equations +% x' = f(x,t), for x(0) = x0 +% over the time interval defined in tspan. +% +% Copyright: +% 2016-2018 - Simo Särkkä and Arno Solin +% +% License: +% This software is provided under the MIT License. See the accompanying +% LICENSE file for details. + +%% + + % Number of steps + steps = numel(tspan); + + % Allocate space + x = zeros(size(x0,1),steps); + + % Initial state + x(:,1) = x0; + + % Iterate + for k=2:steps + + % Time discretization + dt = tspan(k)-tspan(k-1); + + % Stages + dx1 = f(x(:,k-1),tspan(k-1))*dt; + dx2 = f(x(:,k-1)+dx1/2,tspan(k-1)+dt/2)*dt; + dx3 = f(x(:,k-1)+dx2/2,tspan(k-1)+dt/2)*dt; + dx4 = f(x(:,k-1)+dx3,tspan(k-1)+dt)*dt; + + % Step + x(:,k) = x(:,k-1) + 1/6*(dx1 + 2*dx2 + 2*dx3 + dx4); + + end + + \ No newline at end of file diff --git a/matlab/srkS10scalarnoise.m b/matlab/srkS10scalarnoise.m new file mode 100644 index 0000000..5cd34d0 --- /dev/null +++ b/matlab/srkS10scalarnoise.m @@ -0,0 +1,81 @@ +function [x,tspan] = srkS10scalarnoise(f,L,tspan,x0,Q) +%% SRKS10SCALARNOISE - Numerical SDE solver: Stochastic RK, strong 1.0 +% +% Syntax: +% [x,tspan] = srkS10scalarnoise(f,L,tspan,x0,Q) +% +% In: +% f - Drift function, f(x,t) +% L - Diffusion function, L(x,t) +% tspan - Time steps to simulate, [t0,...,tend] +% x0 - Initial condition +% Q - Spectral density (default: standard Brownian motion) +% +% Out: +% x - Solved values +% tspan - Time steps +% +% Description: +% Integrates the system of stochatic differential equations +% dx = f(x,t) dt + L(x,t) dbeta, for x(0) = x0 +% over the time interval defined in tspan. +% +% Copyright: +% 2018 - Simo Särkkä and Arno Solin +% +% License: +% This software is provided under the MIT License. See the accompanying +% LICENSE file for details. + +%% + + % Check if Q given + if nargin<5 || isempty(Q), Q = eye(size(L(x0,tspan(1)),2)); end + + % NB: Only for scalar beta + if size(Q,1)>1, error('NB: Only for scalar beta.'), end + + % Number of steps + steps = numel(tspan); + + % Allocate space + x = zeros(size(x0,1),steps); + + % Initial state + x(:,1) = x0; + + % Pre-calculate random numbers + R = randn(1,steps); + + % Iterate + for k=2:steps + + % Time discretization + dt = tspan(k)-tspan(k-1); + + % Increment + db = sqrt(dt*Q)*R(1,k); + dbb = 1/2*(db^2 - Q*dt); + + % Evaluate only once + fx = f(x(:,k-1),tspan(k-1)); + Lx = L(x(:,k-1),tspan(k-1)); + + % Supporting values + x2 = x(:,k-1) + fx*dt; + tx2 = x2 + Lx*dbb/sqrt(dt); + tx3 = x2 - Lx*dbb/sqrt(dt); + + % Evaluate the remaining values + fx2 = f(x2,tspan(k-1)+dt); + Lx2 = L(tx2,tspan(k-1)+dt); + Lx3 = L(tx3,tspan(k-1)+dt); + + % Step + x(:,k) = x(:,k-1) + ... + (fx+fx2)*dt/2 + ... + Lx*db + ... + sqrt(dt)/2*(Lx2 - Lx3); + + end + diff --git a/matlab/srkW20.m b/matlab/srkW20.m new file mode 100644 index 0000000..84b96ae --- /dev/null +++ b/matlab/srkW20.m @@ -0,0 +1,111 @@ +function [x,tspan] = srkW20(f,L,tspan,x0,Q,gaussian) +%% SRKW20 - Numerical SDE solver: Stochastic RK, weak 2.0 +% +% Syntax: +% [x,tspan] = srkW20(f,L,tspan,x0,Q,gaussian) +% +% In: +% f - Drift function, f(x,t) +% L - Diffusion function, L(x,t) +% tspan - Time steps to simulate, [t0,...,tend] +% x0 - Initial condition +% Q - Spectral density (default: standard Brownian motion) +% gaussian - Use Gaussian increments (default: false) +% +% Out: +% x - Solved values +% tspan - Time steps +% +% Description: +% Integrates the system of stochatic differential equations +% dx = f(x,t) dt + L(x,t) dbeta, for x(0) = x0 +% over the time interval defined in tspan. +% +% Copyright: +% 2018 - Simo Särkkä and Arno Solin +% +% License: +% This software is provided under the MIT License. See the accompanying +% LICENSE file for details. + +%% + + % Number of steps + steps = numel(tspan); + + % Allocate space + x = zeros(size(x0,1),steps); + + % Initial state + x(:,1) = x0; + + % Dimensions + n = size(x0,1); + m = size(L(x0,tspan(1)),2); + + % CHeck inputs + if nargin < 5 || isempty(Q) + Q = eye(m); + end + cQ = chol(Q,'lower'); + + % Pre-calculate random numbers + if nargin < 6 || isempty(gaussian) || ~gaussian + foo = rand(m,steps); + B = sqrt(3)*(foo < 1/6) - sqrt(3)*(foo > 5/6); + foo = rand(m,steps); + Z = -1 + 2*(foo > 1/2); + else + B = randn(m,steps); + foo = rand(m,steps); + Z = -1 + 2*(foo > 1/2); + end + + % Iterate + for k=2:steps + + % Time discretization + dt = tspan(k)-tspan(k-1); + + % Increments + db = cQ*sqrt(dt)*B(:,k); + dbb = (db*db')/2 + ... + sqrt(dt)/2*tril(sqrt(dt)*ones(m,1)*Z(:,k)',-1) - ... + sqrt(dt)/2*triu(sqrt(dt)*Z(:,k)*ones(1,m),1) - ... + dt*eye(m)/2; + + % Evaluate only once + fx = f(x(:,k-1),tspan(k-1)); + Lx = L(x(:,k-1),tspan(k-1)); + + % Supporting values + xx = x(:,k-1) + fx*dt; + tx0 = xx + Lx*db; + Lxdbb = Lx*dbb; + tx2 = bsxfun(@plus,xx,Lxdbb); + tx3 = bsxfun(@plus,xx,-Lxdbb); + + % Iterate to obtain the bar x values + bx2 = tx2; bx3 = tx3; + for n=1:m + ind = setdiff(1:m,n); + bx2(:,n) = xx + sum(Lx(:,ind)*dbb(ind,n),2)/sqrt(dt); %sum(Lxdbb(:,ind),2)/sqrt(dt); + bx3(:,n) = xx - sum(Lx(:,ind)*dbb(ind,n),2)/sqrt(dt); %sum(Lxdbb(:,ind),2)/sqrt(dt); + end + + % Evaluate + Ltx2 = L(tx2,tspan(k)); + Ltx3 = L(tx3,tspan(k)); + Lbx2 = L(bx2,tspan(k)); + Lbx3 = L(bx3,tspan(k)); + + % Make the step + x(:,k) = x(:,k-1) + ... + dt/2*(fx + f(tx0,tspan(k))) + ... + ( Lx/2 + Ltx2/4 + Ltx3/4)*db + ... + (Ltx2/2 - Ltx3/2 )*diag(dbb)/sqrt(dt) + ... + ( -Lx/2 + Lbx2/4 + Lbx3/4)*db + ... + (Lbx2/2 - Lbx3/2 )*ones(m,1)*sqrt(dt); + + end + diff --git a/matlab/w20scalar.m b/matlab/w20scalar.m new file mode 100644 index 0000000..5907993 --- /dev/null +++ b/matlab/w20scalar.m @@ -0,0 +1,97 @@ +function [x,tspan] = w20scalar(f,L,tspan,x0,Q,gaussian) +%% W20scalar - Weak order 2.0 Ito-Taylor +% +% Syntax: +% [x,tspan] = w20scalar(f,L,tspan,x0,Q,gaussian) +% +% In: +% f - Drift function and derivatives, {f, df, ddf} +% L - Diffusion function and derivatives, {L, dL, ddL} +% tspan - Time steps to simulate, [t0,...,tend] +% x0 - Initial condition +% Q - Spectral density (default: standard Brownian motion) +% gaussian - Use Gaussian increments (default: false) +% +% Out: +% x - Solved values +% tspan - Time steps +% +% Description: +% Integrates the system of stochatic differential equations +% dx = f(x,t) dt + L(x,t) dbeta, for x(0) = x0 +% over the time interval defined in tspan. +% +% Copyright: +% 2018 - Simo Särkkä and Arno Solin +% +% License: +% This software is provided under the MIT License. See the accompanying +% LICENSE file for details. + +%% + + % Check if Q given + if nargin<5 || isempty(Q), Q = eye(size(L(x0,tspan(1)),2)); end + + % Check the sort of increments to use + if nargin<6 || isempty(gaussian), gaussian = false; end + + % Extract derivatives + ddf = f{3}; + df = f{2}; + f = f{1}; + ddL = L{3}; + dL = L{2}; + L = L{1}; + + % Cholesky factor of Q + cQ = chol(Q,'lower'); + + % Number of steps + steps = numel(tspan); + + % Allocate space + x = zeros(size(x0,1),steps); + + % Initial state + x(:,1) = x0; + + % Pre-calculate random numbers + RR = rand(size(Q,1),steps); + R = (RR<1/6)*sqrt(3) - (RR>5/6)*sqrt(3); + + % Iterate + for k=2:steps + + % Time discretization + dt = tspan(k)-tspan(k-1); + + % Increment + if gaussian + db = sqrt(dt)*cQ*randn; + dz = 1/2*db*dt; + else + db = sqrt(dt)*cQ*R(:,k); + dz = 1/2*dt^(3/2)*R(:,k); + end + + % Evaluate only once + fx = f(x(:,k-1),[]); + dfx = df(x(:,k-1),[]); + ddfx = ddf(x(:,k-1),[]); + Lx = L(x(:,k-1),[]); + dLx = dL(x(:,k-1),[]); + ddLx = ddL(x(:,k-1),[]); + + % Step + x(:,k) = x(:,k-1) + ... + fx*dt + ... + Lx*db + ... + 1/2*Lx*dLx*(db^2 - dt) + ... + Lx*df(x(:,k-1),[])*dz + ... + 1/2*(fx*dfx + 1/2*Lx^2*ddfx)*dt^2 + ... + (fx*dLx - 1/2*Lx^2*ddLx)*(db*dt - dz); + + end + +