Skip to content

Commit

Permalink
spcr
Browse files Browse the repository at this point in the history
  • Loading branch information
mlesnoff committed Jan 17, 2025
1 parent 6eb7f03 commit ed432c6
Show file tree
Hide file tree
Showing 4 changed files with 15 additions and 26 deletions.
5 changes: 0 additions & 5 deletions src/_structures_fun.jl
Original file line number Diff line number Diff line change
Expand Up @@ -357,14 +357,9 @@ end

struct Spcr
fitm::Spca
T::Matrix
V::Matrix
C::Matrix
xmeans::Vector
xscales::Vector
ymeans::Vector
yscales::Vector
weights::Weight
sellv::Vector{Vector{Int}}
sel::Vector{Int}
par::ParSpcr
Expand Down
10 changes: 4 additions & 6 deletions src/pcr.jl
Original file line number Diff line number Diff line change
Expand Up @@ -80,14 +80,12 @@ function pcr!(X::Matrix, Y::Matrix, weights::Weight; kwargs...)
Q = eltype(X)
q = nco(Y)
ymeans = colmean(Y, weights)
## No need to fscale Y
## below yscales is built only for consistency with coef::Plsr
yscales = ones(Q, q)
## End
yscales = ones(Q, q) # built only for consistency with coef::Plsr
fitm = pcasvd!(X, weights; kwargs...)
## Below, first term of the product is equal to Diagonal(1 ./ fitm.sv[1:nlv].^2)
## if T is D-orthogonal. This is the case for the actual version (pcasvd)
## theta: coeffs regression of Y on T (= C')
## theta: coefs regression of Y on T (= C')
## not needed (same theta): fcenter!(Y, ymeans)
theta = inv(fitm.T' * fweight(fitm.T, fitm.weights.w)) * fitm.T' * fweight(Y, fitm.weights.w)
Pcr(fitm, theta', ymeans, yscales, par)
end
Expand Down Expand Up @@ -118,7 +116,7 @@ function coef(object::Pcr; nlv = nothing)
isnothing(nlv) ? nlv = a : nlv = min(nlv, a)
theta = vcol(object.C, 1:nlv)'
Dy = Diagonal(object.yscales)
## To not use for Spcr (R not computed; while for Pcr, R = V)
## Not used for Spcr (since R not computed; while for Pcr, R = V)
B = fweight(vcol(object.fitm.V, 1:nlv), 1 ./ object.fitm.xscales) * theta * Dy
## In 'int': No correction is needed, since
## ymeans, xmeans and B are in the original scale
Expand Down
24 changes: 10 additions & 14 deletions src/spcr.jl
Original file line number Diff line number Diff line change
Expand Up @@ -53,9 +53,10 @@ nvar = 20
model = spcr(; nlv, meth, nvar) ;
fit!(model, Xtrain, ytrain)
pnames(model)
pnames(model.fitm)
@head model.fitm.T
@head model.fitm.V
fitm = model.fitm ;
pnames(fitm)
@head fitm.fitm.T
@head fitm.fitm.V
@head transf(model, Xtest)
@head transf(model, Xtest; nlv = 3)
Expand Down Expand Up @@ -90,15 +91,10 @@ function spcr!(X::Matrix, Y::Matrix, weights::Weight; kwargs...)
Q = eltype(X)
q = nco(Y)
ymeans = colmean(Y, weights)
## No need to fscale Y
yscales = ones(Q, q)
## End
yscales = ones(Q, q)
fitm = spca!(X, weights; kwargs...)
K = fweight(fitm.T, weights.w)'
beta = inv(K * fitm.T) * K * fcenter(Y, ymeans)
Spcr(fitm, fitm.T, fitm.V, beta', fitm.xmeans, fitm.xscales, ymeans, yscales, weights,
fitm.sellv, fitm.sel, # add compared to ::Pcr
par)
theta = inv(fitm.T' * fweight(fitm.T, fitm.weights.w)) * fitm.T' * fweight(Y, fitm.weights.w)
Spcr(fitm, theta', ymeans, yscales, fitm.sellv, fitm.sel, par)
end

"""
Expand All @@ -110,14 +106,14 @@ Compute Y-predictions from a fitted model.
"""
function predict(object::Spcr, X; nlv = nothing)
X = ensure_mat(X)
a = nco(object.T)
a = nco(object.fitm.T)
isnothing(nlv) ? nlv = a : nlv = (max(0, minimum(nlv)):min(a, maximum(nlv)))
le_nlv = length(nlv)
T = transf(object, X)
pred = list(Matrix{eltype(X)}, le_nlv)
@inbounds for i in eachindex(nlv)
beta = vcol(object.C, 1:nlv[i])'
pred[i] = vcol(T, 1:nlv[i]) * beta .+ object.ymeans'
theta = vcol(object.C, 1:nlv[i])'
pred[i] = vcol(T, 1:nlv[i]) * theta .+ object.ymeans'
end
le_nlv == 1 ? pred = pred[1] : nothing
(pred = pred,)
Expand Down
2 changes: 1 addition & 1 deletion src/vip.jl
Original file line number Diff line number Diff line change
Expand Up @@ -104,7 +104,7 @@ function vip(object::Union{Pcr, Plsr, Spcr}; nlv = nothing)
(imp = imp, W2, sst)
end

function vip(object::Union{Pcr, Plsr}, Y; nlv = nothing)
function vip(object::Union{Pcr, Plsr, Spcr}, Y; nlv = nothing)
if isa(object, Jchemo.Plsr)
W = object.W
T = object.T
Expand Down

0 comments on commit ed432c6

Please sign in to comment.