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.
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.9499009 64.42877
#> True variance of background noise: 20
#> Estimated variance of background noise: 19.84521
#> 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.000272548
#> True variance of background noise: 20
#> Estimated variance of background noise: 20.08427
#> Maximum log likelihood value: -2289361
#> Akaike information criterion score: 4578827
Example 5: Simulate BM and get model-free MSD estimation
set.seed(1)
## Simulation
sim_bm = simulation(sz=100,len_t=100,sigma_bm = 0.5)
sam_mf = SAM(sim_object = sim_bm, model_name = sim_bm@model_name,
model_free_estimation = TRUE)
show(sam_mf)
#> Number of q ring: 49
#> Index of wave number selected: 1 2 3 4 5 6 7 8 9 11 13 15 17 20
#> Number of lag times: 99
#> Maximum log likelihood value: 142743.5
#> MSD uncertainty intervals estimated: FALSE
#> Storage/loss moduli estimated: FALSE
#> True MSD available: TRUE
plot_MSD(object = sam_mf, msd_truth = sam_mf@msd_truth)
Example 6: Real data and get model-free MSD estimation
tif_path = system.file("extdata", "example_video.tif", package = "AIUQ")
## Load the microscopy video which is a space-by-space-by-time array
intensity = bioimagetools::readTIF(tif_path)
## Specify the structure of the input intensity data
intensity_str = "SST_array"
## Provide the time interval between consecutive frames and the pixel size
mindt = 1
pxsz = 1
sam_mf = SAM(intensity = intensity, intensity_str = intensity_str,
mindt = mindt, pxsz = pxsz, model_free_estimation = TRUE)
plot(log10(sam_mf@d_input), log10(sam_mf@msd_est), col = "blue", pch = 16,
xlab = expression(log[10](Delta * t)), ylab = expression(log[10](MSD)))
Example 7: Real data and get storage and loss moduli estimation
## MSD_unit, Temp, and a should be specified according to the experiment.
## Modulus estimation requires uncertainty = TRUE so that MSD intervals can be used
## for sampling. The medians of sampled Gp and Gpp are used as the estimate of
## storage and loss moduli.
sam_mf = SAM(intensity = real_data, intensity_str = "SST_array",
pxsz = pxsz, mindt = mindt,
model_free_estimation = TRUE,
uncertainty = TRUE, M = M,
estimate_moduli = TRUE,
MSD_unit = MSD_unit, Temp = Temp, a = a)
plot(log10(sam_mf@omega), log10(sam_mf@Gp), type = "l", col = "blue",
xlab = expression(log[10](omega)), ylab = expression(log[10](modulus)))
lines(log10(sam_mf@omega), log10(sam_mf@Gpp), col = "red")
legend("topright", legend = c("G'", "G''"), col = c("blue", "red"), lty = 1)