AIUQ tutorial

In this vignette, we will demonstrate the core functionalities of the AIUQ package. These include estimating parameters and mean squared displacement(MSD) with associated uncertainties; simulating particle movements governed by various stochastic processes and generating corresponding intensity profiles to emulate microscopic images. More examples including the application of this package to real experimental data can be found on GitHub.

We start by importing the AIUQ library.

library(AIUQ)

Example 1: Simulate BM and get estimated parameters using BM model

To illustrate the method, we simulate a data set using default values of the simulation class which corresponding to the Brownian Motion(BM). (show() prints the main parameters used in simulation.)

set.seed(1)

sim_bm = simulation()
show(sim_bm)
#> Frame size:  200 200 
#> Number of time steps:  200 
#> Number of particles:  50 
#> Stochastic process:  BM 
#> Variance of background noise:  20 
#> sigma_bm:  1

## Plot simulated particle trajectory
plot_traj(sim_bm)

par(mfrow=c(1,2))
## Plot intensity profile for different frames 
plot_intensity(sim_bm@intensity, sz=sim_bm@sz) #first frame 
plot_intensity(sim_bm@intensity, sz=sim_bm@sz,frame=10, color=T) #10th frame, color image

Next, we can estimate the MSD and other parameters with selected fitting model.

## AIUQ method: use BM as fitted model 
sam = SAM(sim_object=sim_bm, model_name='BM')
show(sam)
#> Fitted model:  BM 
#> Number of q ring:  99 
#> Index of wave number selected:  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 
#> True parameters in the model:  2 
#> Estimated parameters in the model:  2.010371 
#> True variance of background noise:  20 
#> Estimated variance of background noise:  19.90164 
#> Maximum log likelihood value:  -18150726 
#> Akaike information criterion score:  36301654

par(mfrow=c(1,2))
## Plot true MSD and estimated MSD
plot_MSD(object=sam, msd_truth=sam@msd_truth) #in log10 scale

## Plot intensity in reciprocal space 
plot_I_q_1 = matrix(sam@I_q[,1], sam@sz[1],sam@sz[2]) #first frame 
plot3D::image2D(abs(fftshift(plot_I_q_1)),main="intensity in reciprocal space")

User can select wavevector q range via AIUQ_thr or index_q.

sam = SAM(sim_object=sim_bm, AIUQ_thr=c(0.99,0.6)) 
#Note: Default model_name is "BM", it's ok to not specify this argument if want to fit with BM model 
show(sam)
#> Fitted model:  BM 
#> Number of q ring:  99 
#> Index of wave number selected:  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 
#> True parameters in the model:  2 
#> Estimated parameters in the model:  2.008843 
#> True variance of background noise:  20 
#> Estimated variance of background noise:  20.05226 
#> Maximum log likelihood value:  -8065526 
#> Akaike information criterion score:  16131176

sam = SAM(sim_object=sim_bm, index_q_AIUQ=5:50)
show(sam)
#> Fitted model:  BM 
#> Number of q ring:  99 
#> Index of wave number selected:  5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 
#> True parameters in the model:  2 
#> Estimated parameters in the model:  2.015093 
#> True variance of background noise:  20 
#> Estimated variance of background noise:  20.02955 
#> Maximum log likelihood value:  -6198715 
#> Akaike information criterion score:  12397526

Example 2: Simulate BM with user defined parameters and get estimated MSD with uncertainty

set.seed(1)

## Simulation
sim_bm = simulation(sz=100,len_t=100,sigma_bm=0.5)
show(sim_bm)
#> Frame size:  100 100 
#> Number of time steps:  100 
#> Number of particles:  50 
#> Stochastic process:  BM 
#> Variance of background noise:  20 
#> sigma_bm:  0.5

## Plot simulated particle trajectory
plot_traj(sim_bm)

## AIUQ method: fitting using BM model with uncertainty quantification 
sam = SAM(sim_object=sim_bm, uncertainty=T)
show(sam)
#> Fitted model:  BM 
#> Number of q ring:  49 
#> Index of wave number selected:  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 
#> True parameters in the model:  0.5 
#> Estimated parameters in the model:  0.5050874 
#> True variance of background noise:  20 
#> Estimated variance of background noise:  20.08774 
#> Maximum log likelihood value:  -2289361 
#> Akaike information criterion score:  4578825

par(mfrow=c(1,2))
## Plot true MSD and estimated MSD with uncertainty 
plot_MSD(object=sam, msd_truth=sam@msd_truth) #in log10 scale 
plot_MSD(object=sam, msd_truth=sam@msd_truth,log10=F) #in real scale

Example 3: Simulate OU process and get estimated parameters using OU model

set.seed(1)

## Simulation
sim_ou = simulation(sigma_ou=4, model_name="OU")
show(sim_ou)
#> Frame size:  200 200 
#> Number of time steps:  200 
#> Number of particles:  50 
#> Stochastic process:  OU 
#> Variance of background noise:  20 
#> (rho, sigma_ou):  0.95 4

par(mfrow=c(1,2))
## Plot simulated particle trajectory
plot_traj(sim_ou)

## AIUQ method: fitting using OU model
sam_ou = SAM(sim_object=sim_ou, model_name=sim_ou@model_name) 
show(sam_ou)
#> Fitted model:  OU 
#> Number of q ring:  99 
#> Index of wave number selected:  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 
#> True parameters in the model:  0.95 64 
#> Estimated parameters in the model:  0.9499101 64.44048 
#> True variance of background noise:  20 
#> Estimated variance of background noise:  19.84462 
#> Maximum log likelihood value:  -18293754 
#> Akaike information criterion score:  36587711

## Plot true MSD and estimated MSD with uncertainty
plot_MSD(object=sam_ou, msd_truth=sam_ou@msd_truth) #in log10 scale

Example 4: User defined MSD structure

set.seed(1)

## Simulation
sim_bm = simulation(sz=100,len_t=100,sigma_bm=0.5)
show(sim_bm)
#> Frame size:  100 100 
#> Number of time steps:  100 
#> Number of particles:  50 
#> Stochastic process:  BM 
#> Variance of background noise:  20 
#> sigma_bm:  0.5

## User defined MSD structure: function of parameters and 
#                                          vector of lag times 
msd_fn = function(param, d_input){
  MSD = param[1]*d_input+param[2]*d_input^2
}

# show MSD and MSD gradient with a simple example 
theta = c(2,1)
d_input = 0:10
model_name = "user_defined"
MSD_list = get_MSD_with_grad(theta=theta,d_input=d_input,model_name=model_name,
                             msd_fn=msd_fn)
MSD_list$msd
#>  [1]   0   3   8  15  24  35  48  63  80  99 120

## AIUQ method: fitting using user_defined model 
sam = SAM(sim_object=sim_bm, model_name=model_name, msd_fn=msd_fn, num_param=2)
show(sam)
#> Fitted model:  user_defined 
#> Number of q ring:  49 
#> Index of wave number selected:  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 
#> True parameters in the model:  0.5 
#> Estimated parameters in the model:  0.5047796 0.0002725495 
#> True variance of background noise:  20 
#> Estimated variance of background noise:  20.08427 
#> Maximum log likelihood value:  -2289361 
#> Akaike information criterion score:  4578827