Skip to content

Commit

Permalink
First take on pacf
Browse files Browse the repository at this point in the history
  • Loading branch information
andreasnoack committed Oct 12, 2013
1 parent c3ec659 commit 72e52c6
Show file tree
Hide file tree
Showing 5 changed files with 136 additions and 3 deletions.
1 change: 1 addition & 0 deletions README.md
Original file line number Diff line number Diff line change
Expand Up @@ -22,6 +22,7 @@ Basic statistics functions for Julia
* `midrange(a)`: Compute the mid point of the range of `a` (e.g `(max(a) + min(a) / 2)`).
* `modes(a)`: Compute all modes of `a`. Be warned that every element of an array with no repeated elements is considered a mode.
* `indicators(a)`: Encode categories using one-hot scheme aka one-of-C encoding, indicator matrix or dummy variables. Optionally, you can provide a list of possible values, e.g. ["A", "B, "C"] or [1:3].
* `pacf(A, l, method)`: Compute partial autocorrelation of an array `A` at lag(s) `l`. The computational method acn either be :regression (default) or :yulewalker.
* `percentile(a)`: Compute the percentiles (0%, 10%, ..., 100%) of `a`.
* `quantile(a)`: Compute any desired quantile of `a`.
* `quartile(a): Compute the quartiles of `a`.
Expand Down
8 changes: 6 additions & 2 deletions src/Stats.jl
Original file line number Diff line number Diff line change
Expand Up @@ -41,21 +41,24 @@ module Stats
acf,
ccf,
cor_spearman,
cor_kendall,
cor_kendall,
pacf,

# others
## Types
StatisticalModel,
RegressionModel,

## Fcuntions
## Functions
coef,
coeftable,
confint,
durbin,
ecdf,
findat,
indicators,
inverse_rle,
levinson,
loglikelihood,
nobs,
predict,
Expand All @@ -68,6 +71,7 @@ module Stats
include("means.jl")
include("scalar_stats.jl")
include("intstats.jl")
include("toeplitzsolvers.jl")
include("corr.jl")
include("others.jl")

Expand Down
48 changes: 47 additions & 1 deletion src/corr.jl
Original file line number Diff line number Diff line change
Expand Up @@ -234,4 +234,50 @@ ccf{T<:Real}(X::AbstractMatrix{T}, args1...; args2...) = ccf(float(X), args1...;
ccf{T<:Real}(x::AbstractMatrix{T}, lags::Integer; correlation::Bool=true, demean::Bool=true) =
reshape(ccf(x, lags:lags, correlation=correlation, demean=demean), size(x, 2), size(x, 2))


# Partial autoroccelation
function pacf{T<:BlasReal}(X::AbstractMatrix{T}, lags::AbstractVector{Int} = 0:min(size(X,1)-1, int(10log10(size(X,1))));
method::Symbol = :regression)
n, p = size(X)
nk = length(lags)
mk = max(lags)
if min(lags) < 0 error("Negative autoroccelations not allowed") end
if 2mk >= n error("Can at most calculate pacf for $(div(n,2) - 1) lags, you requested $mk") end
val = Array(T, nk, p)
if method == :regression
tmpX = ones(T, n, mk + 1)
for j = 1:p
for l = 1:mk
for i = 1+l:n
tmpX[i,l+1] = X[i-l,j]
end
end
i = 1
for l in lags
sX = sub(tmpX, 1+l:n, 1:l+1)
val[i,j] = (cholfact!(sX'sX)\(sX'sub(X, 1+l:n, j)))[end]
i += 1
end
end
elseif method == :yulewalker
tmp = Array(T, mk)
for j = 1:p
acfs = acf(sub(X,1:n,j),1:mk)
i = 1
for l in lags
if l == 0
val[i,j] = one(T)
elseif l == 1
val[i,j] = acfs[i]
else
val[i,j] = -durbin!(sub(acfs, 1:l), tmp)[l]
end
i += 1
end
end
else
error("No such method")
end
return val
end
pacf{T<:Real}(X::AbstractMatrix{T}, args1...; args2...) = pacf(float(X), args1...; args2...)
pacf(x::AbstractVector, args1...; args2...) = squeeze(pacf(reshape(x, length(x), 1), args1...; args2...),2)
66 changes: 66 additions & 0 deletions src/toeplitzsolvers.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,66 @@
# Symmetric Toeplitz solver
function durbin!{T<:BlasReal}(r::AbstractVector{T}, y::AbstractVector{T})
n = length(r)
n <= length(y) || throw(DimensionMismatch("Auxiliary vector cannot be shorter than data vector"))
y[1] = -r[1]
β = one(T)
α = -r[1]
for k = 1:n-1
β *= one(T) - α*α
α = -r[k+1]
for j = 1:k
α -= r[k-j+1]*y[j]
end
α /= β
for j = 1:div(k,2)
tmp = y[j]
y[j] += α*y[k-j+1]
y[k-j+1] += α*tmp
end
if isodd(k) y[div(k,2)+1] *= one(T) + α end
y[k+1] = α
end
return y
end
durbin{T<:BlasReal}(r::AbstractVector{T}) = durbin!(r, zeros(T, length(r)))

function levinson!{T<:BlasReal}(r::AbstractVector{T}, b::AbstractVector{T}, x::AbstractVector{T})
n = length(b)
n == length(r) || throw(DimensionMismatch("Vectors must have same length"))
n <= length(x) || throw(DimensionMismatch("Auxiliary vector cannot be shorter than data vector"))
x[1] = b[1]
b[1] = -r[2]/r[1]
β = one(T)
α = -r[2]/r[1]
for k = 1:n-1
β *= one(T) - α*α
μ = b[k+1]
for j = 2:k+1
μ -= r[j]/r[1]*x[k-j+2]
end
μ /= β
for j = 1:k
x[j] += μ*b[k-j+1]
end
x[k+1] = μ
if k < n - 1
α = -r[k+2]
for j = 2:k+1
α -= r[j]*b[k-j+2]
end
α /= β*r[1]
for j = 1:div(k,2)
tmp = b[j]
b[j] += α*b[k-j+1]
b[k-j+1] += α*tmp
end
if isodd(k) b[div(k,2)+1] *= one(T) + α end
b[k+1] = α
end
end
for i = 1:n
x[i] /= r[1]
end
return x
end
levinson{T<:BlasReal}(r::AbstractVector{T}, b::AbstractVector{T}) = levinson!(r, copy(b), zeros(T, length(b)))
16 changes: 16 additions & 0 deletions test/01.jl
Original file line number Diff line number Diff line change
Expand Up @@ -180,6 +180,22 @@ for i =1:length(racvx12nodemean)
@test_approx_eq racvx12nodemean[i] jacvx12nodemean[i]
end
println("OK")
print("Test partial autocorrelation (regression): ")
@test_approx_eq pacf(x[:,1], 1:4) [-0.218158122381419,
0.195015316828711,
0.144315804606139,
-0.199791229449779]



println("OK")

print("Test partial autocorrelation (Yule Walker): ")
@test_approx_eq pacf(x[:,1], 1:4, method = :yulewalker) [-0.221173011668873,
0.189683314308021,
0.111857020733719,
-0.175020669835420]
println("OK")

@test iqr([1, 2, 3, 4, 5]) == [2.0, 4.0]

Expand Down

0 comments on commit 72e52c6

Please sign in to comment.