| Title: | Software to Perform Haar Fisz Transforms |
|---|---|
| Description: | A Haar-Fisz algorithm for Poisson intensity estimation. Will denoise Poisson distributed sequences where underlying intensity is not constant. Uses the multiscale variance-stabilization method called the Haar-Fisz transform. Contains functions to carry out the forward and inverse Haar-Fisz transform and denoising on near-Gaussian sequences. Can also carry out cycle-spinning. Main reference: Fryzlewicz, P. and Nason, G.P. (2004) "A Haar-Fisz algorithm for Poisson intensity estimation." Journal of Computational and Graphical Statistics, 13, 621-638. <doi:10.1198/106186004X2697>. |
| Authors: | Piotr Fryzlewicz [aut], Guy Nason [cre] |
| Maintainer: | Guy Nason <[email protected]> |
| License: | GPL (>= 2) |
| Version: | 4.5.4 |
| Built: | 2024-12-24 06:34:53 UTC |
| Source: | CRAN |
Package to denoise Poisson distributed sequence where underlying intensity is not constant. Uses the multiscale variance-stabilization method called the Haar-Fisz transform. Contains functions to carry out the foward and inverse Haar-Fisz transform and denoising on near-Gaussian sequences. Can also carry out cycle-spinning.
Package to denoise Poisson distributed sequence where
underlying intensity is not constant. Uses the multiscale
variance-stabilization method called the Haar-Fisz transform.
Contains functions to carry out the foward and inverse
Haar-Fisz transform and denoising on near-Gaussian sequences.
Can also carry out cycle-spinning.
See main routine denoise.poisson
Piotr Fryzlewicz>
Fryzlewicz, P. (2003) Wavelet Techniques for Time Series and Poisson Data. PhD Thesis, University of Bristol, Bristol, UK https://www.ma.imperial.ac.uk/~gnason/Research/MAPZFthesis.ps.gz
Fryzlewicz, P. and Nason, G.P. (2004) A Haar-Fisz algorithm for Poisson intensity estimation. Journal of Computational and Graphical Statistics, 13, 621-638. doi:10.1198/106186004X2697
Nason, G.P. (2008) Wavelet Methods in Statistics with R. Springer: New York (Section 6.4) doi:10.1007/978-0-387-75961-6
# # # Main Poisson denoising function is denoise.poisson # # Forward Haar-Fisz transform is hft # # Inverse Haar-Fisz transform is hft.inv# # # Main Poisson denoising function is denoise.poisson # # Forward Haar-Fisz transform is hft # # Inverse Haar-Fisz transform is hft.inv
Main routine of the package. Estimates the deterministic discretised intensity of a one-dimensional Poisson process using the Haar-Fisz transformation and partial cycle spinning.
denoise.poisson(y, meth.1 = hf.bt, cs.1 = 50, meth.2 = hf.cv, cs.2 = 50, hybrid = TRUE)denoise.poisson(y, meth.1 = hf.bt, cs.1 = 50, meth.2 = hf.cv, cs.2 = 50, hybrid = TRUE)
y |
The vector of Poisson counts, its length must be a power of 2. |
meth.1 |
Unquoted name of an S-Plus routine for denoising Gaussian contaminated vectors. Must take and return a vector of length |
cs.1 |
The number of cycle spins to be performed with |
meth.2 |
Of the same type as |
cs.2 |
The number of cycle spins to be performed with |
hybrid |
If set to TRUE, then the estimates are computed using both |
For a given input sequence, basic operation of the code
performs a cyclic shift on
the data. Then applies the Haar-Fisz transform, then one
of the denoising methods (specified by meth.1),
then the inverse Haar-Fisz transform and then a shift back.
This is repeated for cs.1 cyclic shifts and the results of
all shifts returned.
Returns vector of the same length as the input y but is the
denoised estimate.
Piotr Fryzlewicz
Fryzlewicz, P. and Nason, G.P. (2004) A Haar-Fisz algorithm for Poisson intensity estimation. Journal of Computational and Graphical Statistics, 13, 621-638. doi:10.1198/106186004X2697
hft, hft.inv,
hf.u, hf.cv, hf.bt, hf.tiu
# # Apply denoise.poisson to xquake data # data(xquake) xquake.denoised <- denoise.poisson(xquake) # # Now plot the original data and it's denoised version in red # plot(xquake, type="l") lines(xquake.denoised, col=2)# # Apply denoise.poisson to xquake data # data(xquake) xquake.denoised <- denoise.poisson(xquake) # # Now plot the original data and it's denoised version in red # plot(xquake, type="l") lines(xquake.denoised, col=2)
Denoises a Gaussian contaminated vector using a version of the wavelet-based “greedy tree" algorithm by Baraniuk, see references. Note: this function does not actually do any Haar-Fisz transform, it is a homogeneous variance Gaussian denoising function, which is used by the haarfisz package.
hf.bt(x, filter.number = 1, family = "DaubExPhase", min.level = 3, noise.level = NULL)hf.bt(x, filter.number = 1, family = "DaubExPhase", min.level = 3, noise.level = NULL)
x |
The noisy vector, its length must be a power of 2. |
filter.number |
The filter number of the analysing wavelet. Can be set to 1, 2, ..., 10 for |
family |
The family of wavelet bases from which the wavelet |
min.level |
The minimum level thresholded. |
noise.level |
Standard deviation of the noise, can be set to a positive number or NULL; in the latter case it will be estimated using MAD. |
Denoised version of x.
Piotr Fryzlewicz
Baraniuk, R. G. (1999) Optimal tree approximation with wavelets. Pages 206-214 of Unser, M.A., Aldroubi, A. and Laine, A.F. (eds), Wavelet Applications in Signal and Image Processing VII, Proceedings of SPIE 3813. SPIE. doi:10.1117/12.366780
Fryzlewicz, P. and Nason, G.P. (2004) A Haar-Fisz algorithm for Poisson intensity estimation. Journal of Computational and Graphical Statistics, 13, 621-638. doi:10.1198/106186004X2697
# # Generate simple sinusoidal test signal # test.sig <- sin(seq(from=0, to=6*pi, length=128)) # # Invent simulated noisy signal # test.dat <- test.sig + rnorm(128, sd=0.2) # # Denoise using hf.bt # test.est <- hf.bt(test.dat) # # Now plot the results: the truth, the noisy signal, the estimate # ts.plot(test.dat) lines(test.est, col=2) lines(test.sig, col=3, lty=2)# # Generate simple sinusoidal test signal # test.sig <- sin(seq(from=0, to=6*pi, length=128)) # # Invent simulated noisy signal # test.dat <- test.sig + rnorm(128, sd=0.2) # # Denoise using hf.bt # test.est <- hf.bt(test.dat) # # Now plot the results: the truth, the noisy signal, the estimate # ts.plot(test.dat) lines(test.est, col=2) lines(test.sig, col=3, lty=2)
Denoises a Gaussian contaminated vector using wavelet thresholding with a threshold chosen by “leave-half-out" cross-validation. Note: this function does not actually do any Haar-Fisz transform, it is a homogeneous variance Gaussian denoising function, which is used by the haarfisz package.
hf.cv(x, filter.number = 1, family = "DaubExPhase", min.level = 3, type = "hard")hf.cv(x, filter.number = 1, family = "DaubExPhase", min.level = 3, type = "hard")
x |
The noisy vector, its length must be a power of 2. |
filter.number |
The filter number of the analysing wavelet. Can be set to 1, 2, ..., 10 for |
family |
The family of wavelet bases from which the wavelet |
min.level |
The minimum level thresholded. |
type |
Type of thresholding, can be set to "hard" or "soft". |
Uses threshold,
wd and
AvBasis
Denoised version of x
Piotr Fryzlewicz
Fryzlewicz, P. and Nason, G.P. (2004) A Haar-Fisz algorithm for Poisson intensity estimation. Journal of Computational and Graphical Statistics, 13, 621-638. doi:10.1198/106186004X2697
# # Generate simple sinusoidal test signal # test.sig <- sin(seq(from=0, to=6*pi, length=128)) # # Invent simulated noisy signal # test.dat <- test.sig + rnorm(128, sd=0.2) # # Denoise using hf.bt # test.est <- hf.cv(test.dat) # # Now plot the results: the truth, the noisy signal, the estimate # ts.plot(test.dat) lines(test.est, col=2) lines(test.sig, col=3, lty=2)# # Generate simple sinusoidal test signal # test.sig <- sin(seq(from=0, to=6*pi, length=128)) # # Invent simulated noisy signal # test.dat <- test.sig + rnorm(128, sd=0.2) # # Denoise using hf.bt # test.est <- hf.cv(test.dat) # # Now plot the results: the truth, the noisy signal, the estimate # ts.plot(test.dat) lines(test.est, col=2) lines(test.sig, col=3, lty=2)
Denoises a Gaussian contaminated vector using translation-invariant hard wavelet thresholding with the universal threshold. Note: this function does not actually do any Haar-Fisz transform, it is a homogeneous variance Gaussian denoising function, which is used by the haarfisz package.
hf.tiu(x, filter.number = 1, family = "DaubExPhase", min.level = 3, noise.level = 1)hf.tiu(x, filter.number = 1, family = "DaubExPhase", min.level = 3, noise.level = 1)
x |
The noisy vector, its length must be a power of 2. |
filter.number |
The filter number of the analysing wavelet. Can be set to 1, 2, ..., 10 for |
family |
The family of wavelet bases from which the wavelet |
min.level |
The minimum level thresholded. |
noise.level |
Standard deviation of the noise, can be set to a positive number or to an estimate (a function of x). |
Uses threshold,
wd and
AvBasis
Denoised version of x.
Piotr Fryzlewicz
Fryzlewicz, P. and Nason, G.P. (2004) A Haar-Fisz algorithm for Poisson intensity estimation. Journal of Computational and Graphical Statistics, 13, 621-638. doi:10.1198/106186004X2697
# # Generate simple sinusoidal test signal # test.sig <- sin(seq(from=0, to=6*pi, length=128)) # # Invent simulated noisy signal # test.dat <- test.sig + rnorm(128, sd=0.2) # # Denoise using hf.bt # test.est <- hf.tiu(test.dat) # # Now plot the results: the truth, the noisy signal, the estimate # ts.plot(test.dat) lines(test.est, col=2) lines(test.sig, col=3, lty=2)# # Generate simple sinusoidal test signal # test.sig <- sin(seq(from=0, to=6*pi, length=128)) # # Invent simulated noisy signal # test.dat <- test.sig + rnorm(128, sd=0.2) # # Denoise using hf.bt # test.est <- hf.tiu(test.dat) # # Now plot the results: the truth, the noisy signal, the estimate # ts.plot(test.dat) lines(test.est, col=2) lines(test.sig, col=3, lty=2)
Denoises a Gaussian contaminated vector using wavelet thresholding with the universal threshold. Note: this function does not actually do any Haar-Fisz transform, it is a homogeneous variance Gaussian denoising function, which is used by the haarfisz package.
hf.u(x, filter.number = 10, family = "DaubLeAsymm", min.level = 3, type = "hard")hf.u(x, filter.number = 10, family = "DaubLeAsymm", min.level = 3, type = "hard")
x |
The noisy vector, its length must be a power of 2. |
filter.number |
The filter number of the analysing wavelet. Can be set to 1, 2, ..., 10 for |
family |
The family of wavelet bases from which the wavelet |
min.level |
The minimum level thresholded. |
type |
Type of thresholding, can be set to hard or soft |
Uses threshold,
wd and
AvBasis
Denoised version of x.
Piotr Fryzlewicz
Fryzlewicz, P. and Nason, G.P. (2004) A Haar-Fisz algorithm for Poisson intensity estimation. Journal of Computational and Graphical Statistics, 13, 621-638. doi:10.1198/106186004X2697
# # Generate simple sinusoidal test signal # test.sig <- sin(seq(from=0, to=6*pi, length=128)) # # Invent simulated noisy signal # test.dat <- test.sig + rnorm(128, sd=0.2) # # Denoise using hf.bt # test.est <- hf.u(test.dat) # # Now plot the results: the truth, the noisy signal, the estimate # ts.plot(test.dat) lines(test.est, col=2) lines(test.sig, col=3, lty=2)# # Generate simple sinusoidal test signal # test.sig <- sin(seq(from=0, to=6*pi, length=128)) # # Invent simulated noisy signal # test.dat <- test.sig + rnorm(128, sd=0.2) # # Denoise using hf.bt # test.est <- hf.u(test.dat) # # Now plot the results: the truth, the noisy signal, the estimate # ts.plot(test.dat) lines(test.est, col=2) lines(test.sig, col=3, lty=2)
Performs the (forward) Haar-Fisz transform.
hft(data)hft(data)
data |
A vector of Poisson counts, its length must be a power of 2 |
The Haar-Fisz for Poisson works, roughly speaking, by taking
the Haar wavelet transform of data. Then dividing the
mother wavelet coefficients by the respective father coefficients,
and replacing the results of the divisions back into the same
coefficient locations, and then carrying out an inverse
Haar wavelet transform. This produces a nearer-Gaussian variance
stabilized version of the original (or a version of the underlying
intensity which is close to an 'intensity PLUS homogeneous
Gaussian noise' signal, which is easier to denoise using ‘standard’
methods.
The inverse transform is hft.inv
The Haar-Fisz transform of data, which will be the same
length as data.
Piotr Fryzlewicz
Fryzlewicz, P. and Nason, G.P. (2004) A Haar-Fisz algorithm for Poisson intensity estimation. Journal of Computational and Graphical Statistics, 13, 621-638. doi:10.1198/106186004X2697
denoise.poisson,
hft.inv,
hf.bt,
hf.cv,
hf.u,
hf.tiu
# # Generate Poisson data, half with one intensity, and half with a larger one # v <- c( rpois(64, lambda=1), rpois(64, lambda=10)) # # Plot it to note that the variation is bigger in the second half # (and the mean, but this is not important for this bit) # ts.plot(v) # # Now do the Haar-Fisz transform # vhft <- hft(v) # # Now plot this, and see that the variance of the second bit is now comparable # to the first # ts.plot(vhft)# # Generate Poisson data, half with one intensity, and half with a larger one # v <- c( rpois(64, lambda=1), rpois(64, lambda=10)) # # Plot it to note that the variation is bigger in the second half # (and the mean, but this is not important for this bit) # ts.plot(v) # # Now do the Haar-Fisz transform # vhft <- hft(v) # # Now plot this, and see that the variance of the second bit is now comparable # to the first # ts.plot(vhft)
Performs the inverse Haar-Fisz transform.
hft.inv(data)hft.inv(data)
data |
Vector of length |
The inverse Haar-Fisz transform of data
(vector of the same length as data).
Piotr Fryzlewicz
Fryzlewicz, P. and Nason, G.P. (2004) A Haar-Fisz algorithm for Poisson intensity estimation. Journal of Computational and Graphical Statistics, 13, 621-638. doi:10.1198/106186004X2697
# # Make up test set (mimics sequence with half low intensity, followed by # half high intensity) # test.set <- c(8,5,6,3, 30,40,20,35) # # Do forward HFT # test.hft <- hft(test.set) test.hft # [1] 16.38621 15.20951 15.65216 14.23795 21.20421 22.89452 19.27753 22.13791 # # Do inverse HFT # test.back <- hft.inv(test.hft) test.back # [1] 8 5 6 3 30 40 20 35 # # Same as original ## # Make up test set (mimics sequence with half low intensity, followed by # half high intensity) # test.set <- c(8,5,6,3, 30,40,20,35) # # Do forward HFT # test.hft <- hft(test.set) test.hft # [1] 16.38621 15.20951 15.65216 14.23795 21.20421 22.89452 19.27753 22.13791 # # Do inverse HFT # test.back <- hft.inv(test.hft) test.back # [1] 8 5 6 3 30 40 20 35 # # Same as original #
One of my functions to resolve issues for a similar function that seems to have been forgotten in haarfisz.
shift.sequence(v, places, dir="right")shift.sequence(v, places, dir="right")
v |
Vector to shift |
places |
The number of places to shift |
dir |
Whether the shift should be right or left |
This function takes a sequence input and shifts it to the left or right by the specified number of places. This is a circular shift. For example, when shifting to the right, any numbers that drop off are appended circularly to the front, etc.
a shifted output sequence.
Piotr Fryzlewicz
Fryzlewicz, P. and Nason, G.P. (2004) A Haar-Fisz algorithm for Poisson intensity estimation. Journal of Computational and Graphical Statistics, 13, 621-638. doi:10.1198/106186004X2697
# # Shift 1:10 one place to the right # shift.sequence(1:10,1, dir="right") # [1] 10 1 2 3 4 5 6 7 8 9 # # Shift 1:10 twos place to the right # shift.sequence(1:10,2, dir="right") # [1] 9 10 1 2 3 4 5 6 7 8 # # Shift 1:10 one place to the left # shift.sequence(1:10,1, dir="left") # [1] 2 3 4 5 6 7 8 9 10 1# # Shift 1:10 one place to the right # shift.sequence(1:10,1, dir="right") # [1] 10 1 2 3 4 5 6 7 8 9 # # Shift 1:10 twos place to the right # shift.sequence(1:10,2, dir="right") # [1] 9 10 1 2 3 4 5 6 7 8 # # Shift 1:10 one place to the left # shift.sequence(1:10,1, dir="left") # [1] 2 3 4 5 6 7 8 9 10 1
Time series of the number of earthquakes of magnitude >= 3.0 which occurred in Northern California in 1024 weeks, the last week being 29/11 – 05/12/2000.
data(xquake)data(xquake)
The series was composed using the data obtained from the Northern California Earthquake Data Center.